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ABSTRACT 


Four methods for determining the geometric elements 
Olwat eelipsing Dbitary from Je dignt curve are vexplored 
in detail. The methods discussed are those of Russell 
(specifically, the version due to Tabachnik), Kitamura, 
Kopal (frequency domain approach), and Wood. In each case, 
the underlying model of an eclipsing binary system is dis- 
cussed. The various methods of light analysis are then 
applied to the eclipsing binaries HS Herculis, W Delphini, 
and HD 219634. The results of each analysis are discussed, 
and the various methods of analysis are compared with one 
another. Finally, the relative merits of each model of an 
eclipsing binary system are considered. Computer programs 
for the various methods of light curve analysis, along with 


explanations of their use, are presented in the appendices. 
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CHAPTER aL 
INTRODUCTION 


Eclipsing binary stars are, in many ways, very 
informative to the astronomer and astrophysicist. The 
study of eclipsing binary stars can reveal a great deal 
abouts the sizes Of Stars, their masses, densities, and 
PiLerials Structure. <SCuch information 1S also: valuable an 
ascertaining the evolutionary state of the two (or more) 
Stal so) CONStltubLing an eclipsing binary system. —“ilm certain 
cases in which the two component stars of an eclipsing 
binary system are in close proximity to one another, it is 
also possible to deduce the amount of tidal distortion 
present, and to check for the presence of matter streams 
between the two stars. With all of this information in 
hand, a scale model of an eclipsing binary star may be 
constructed, and hence, our knowledge of the system will 
be complete. One may then use this scale model to look for 
long-term effects such as apsidal motion, the presence of 
which can be deduced from a accnilar Saeaee in the period of 
the eclipsing binary. Such measurements may also be used 
to verify Einstein's theory of General Relativity. 

The key problem, however, is in interpreting the 
observed light changes (the "light curve') of an eclipsing 
Bbinaby. ) lite 1sethe criticalesecep whilel -kves betweensmalk— 
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in the previous paragraph. To this end a great deal of 
work has been done, from the first tentative steps taken 
by Russelbeingl912 WRusseli,) LOV2)ssewhich dealt withthe 
determination of the geometric elements (the relative radii 
CieeheGstars,, thevorbrtal sinelinagtion angle p~sand™ eer re la 
tive luminosttses) of van teclipsimd binary "consisting of 
non-limb darkened spherical stars, to the recent Fourier 
analysis methods of Kopal (see for instance, Kopal, 1979), 
which use the harmonic content of the observed light 
changes to deduce the geometric elements. 

A necessary ingredient in all methods for determin- 
ing the geometric elements from the observed light changes 
is a realistic physical model of the binary star, preferably 
involving as few assumptions as possible regarding the 
forms of the stars (i.e., spherical, non-spherical) and 
their physical properties (temperature, luminosity, etc.). 
The complexity, and consequently the realism, of such models 
of eclipsing binaries has grown since the initial investi- 
Gavron Of ehe sprobwiemiby sRusseiis ino] 2s Poesent iday 
models of “eclipsing: binaries descrabe a range Of ‘situations, 
from that of two well-separated spherical stars to systems 
ime wich bOciwestabs ace in COntece, im watch case beta Stars 
areweneatly distornred by Eheim micual tidal einteraceton. 

Be question Of some importance is) Chen: which ot the several 
models on eclipsing binery Stars ecuULEecntly  wavad lable bese 
describes a Given eclipsing binary stare To answer this 


question, representative models of eclipsing binary stars 


neéet eaete 6 giteines get + 
pelea seek 
=slex ont mis ee paneer. 
t0 ndttalends yaanie Heiestis ioc 


geranor sieges ers <5 ergh At tos 
abe 


(rer ,beqat-.osustenl ee cS 
dipkt Fevaaedis a 4b ohdacs oni td : 
fesuemnls or: Sontag are engbsh G2 at ao 
epinastes: tes ebodsometic hl scbitot vussessan & a 
aspséss Jfpih beaseqo aa) mova atittarsio 3LTIaOee oAt pak. - 
Yfdexs3ei1g “late Yishid ett sa Tstom £SSkaNiG oise ste = a 
eri oni ys wa S{dieeod a Hnlrahurrs wet 26 valved 
nf een steers "3 ) -exade wits ig sarre’ 
Saplexeans) asad ganate et soaetr ziada 


bas. (£enizete 
ede .yoteoacme hy 
uisten fews 19 iderlect.* 17 vi sqaet genet: hed ov eeiqecs ert 
—cdesitis [bis ink. Sds denise wor, tcc ebiisn so seaeqgitos to 
Yeh, Sageert ‘9 STO rs irbsent Fr pisdacttey! eRe” Be notte 
eupiesutie to epus: » otiavesb. 2eiteaig wide tm +c eliebom 
emadewe oF Giste Leo sscige ceiserxssiee=f ows. 20 toch? seo? 
axeta fad sezoutpitte ni \fastoos miposs eaace Atod firetled ) ok a 
.@otwvtetevnl’ dene becoim tiers a bassoseth ylseo36 wa > 
tdsaves Sa2 te/dntaw het ai) sodeseeamt emog to onkasdip Ae 
gade ofdellava Visasiows- 21892 yte02d. enzegeios Xo. ———_ : 
elds t6vans of SrESe edit: preva hing: nevig \€ 


erass yxsctc guivasiss Yo afehdn avi streets ate 


and their accompanying methods of light curve analysis will 
have to be analyzed and the appropriate conclusions drawn. 
There are three broad categories of methods used in 
light curve analysis. They are the "classical" or "Russell- 
type" methods previously referred to, the "synthesis" 
methods, devised in the early 1970s, and the "frequency- 
domain" methods of Kopal previously referred to. Two 
methods of the "classical" type are the Russell-Merrill 
(1952) method and the method of Kitamura. A fine example 
of a "synthesis" method is a method devised by Wood. Seve- 
ral versions of Kopal's frequency-domain method exist, but 
the best of these are the most recent ones (e.g., Kopal 
(1982)). These methods of light curve analysis will be 
applied to three stars covering a wide range of physical 
conditions, from the well-separated case of W Delphini, 
to Ene case Of HS »pHerculis with its matter stream, and 
finally to HD 219634, which may be a massive binary and 
possibly even an X-ray source (see Gulliver, Hube and Lowe 
(1982)). This analysis will, we hope, shed some light on 
the validaty Of “the various models Jor eclipsing binery Ssvars 


to be considered. 
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CHAPTER: 2 
THE RUSSELL MODEL 


2.1 General Principles 

The first steps toward an interpretation of eclips- 
ing binary light curves were taken by Russell in 1912 
(Russell (1912a,b)). Subsequent refinements to the theory 
were made by Russell and two collaborators, Merrill and 
Shapley (see Russell and Shapley (1912); Russell and Merrill 
(1952)). The Russell model can be applied to both spherical 
and non-spherical stars with varying degrees of accuracy. 

The TSRINECA CES, aS assumes that the eclipsing 
binary system consists of two spherical stars moving in 
ClEcuvar OEDLUS about a Common Centre Of gravity. The dis— 
tribution of surface brightness J(y) over the disk of each 


star is assumed to follow the "cosine law" 
WEG) =a 0) = eee Se eCOs aya) G2 18) 


where J(0) is the surface brightness at the centre of the 
observed disk of either star, x the coefficient of limb 
darkening, and y the angle of foreshortening, or the angle 
between the line of sight and a radius vector from the 
centre of the star (see figure 1). The angle y varies 
between 0 and 90 degrees. At this point, it should be 
Noted that tne physical properties ob the starseenter tne 
Russell model only through equation, (2.1). The detailed 


features Ob the Stars (e.9., Starspote, Magnetic ticids, 
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FIGURE 1, DEFINITION OF Y, THE ANGLE OF 
FORESHORTENING, 
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TO arrive at a method for determining the geometric 
elements of an eclipsing binary using the spherical model, 
it will be necessary to consider, in some detail, the 
geometry of an eclipse. The treatment that follows is 
Similars Go that given by Irwin (1962) or Russell (191Za,b).. 
Before plunging headlong into this problem, it will be 
necessary to define some of the quantities that will be 
used. The radius of the smaller star is denoted by ter and 
that. Of they larger star Dy ro- These radii are measured in 
units of the centre-to-centre separation of the component 
stars of an eclipsing binary. The relative luminosity of 
the small and large stars will be denoted by Le and fe 


respectively. These luminosities are so defined that 
ipo C= eee: UZe2) 


Since an eclipsing binary light curve displays brightness 
as a function of time, it is necessary to define an orbital 


phase 9 by 


where P is the period of the eclipsing binary, t the time 

at which the brightness was observed or is to be calculated, 
and ty the time of conjunction, which usually coincides with 
the time of minimum light during the primary (deeper) 
eclipse. The times t and oe are measured in Julian days, 


while the period P is measured in days. With explicit 
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retrerence to the eclipse, “le is customary to detine three 
additional quantities, namely k, p, and 6. The dimension-— 
fess parameter kK is Simply the ratio of the radia re and 


es 
| 


The "geometric depth" p is another dimensionless parameter, 
which represents the extent to which the eclipse has pro- 
gressed at any eclipse phase 98. A parameter closely 
related to p is $6, the apparent separation of the centres 
of the disks of the stars. The parameters p and 6 are 
shown in figure 2. The following equation gives the rela- 


tionship between p and 6 (and does, in fact, serve to 


define p) 
Oname 
jo = songs! (2532) 
Ss 
pe 
§ = et Son) ne (2h si) 


using the definition of k stated above. Relation (2.3b) is 
the more useful of the two relations relating p and ¢; The 
quantity 6 may also be related to the phase @ and the orbi- 
tal inclination. The Orbical inclineticone1s detlined as tne 
angle between a plane perpendicular to the line of sight 
(the “celestial sphere") and the orbital plane (see figure 
3). Through the use of some simple trigonometry, and 
Pecallingvenat the EwoOustans COnsti cueing senewecliosing 
Dinlary Nave a Und Separation,, one arrives aim tne, -geometri— 
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Equation (2.4) is the fundamental equation of all Russell- 


type methods for the determination of the geometric elements. 


The derivation of the geometric relation is outlined in 
figure 4. “It should™@also be noted that only a relative 
orbit is considered, namely, the relative orbit of the 
smaller star about the larger one. In the case of an ellip- 
tical “orbit, the geometric relation would be multiplied by 
ae R being the separation between the stars at any orbital 
phase. The orbital phase would have to be replaced by v-4, 
where v is the true anomaly and w the angle between the 
lamer oresapsides and sihe line cl sight. —ihe vase majoricy 
of eclipsing binary systems have circular orbits, however, 
largely as a consequence of their short orbital periods and 
consequent tidal interations. 

As an aid in the interpretation of eclipsing binary 
Migne curves, a relative luminosity ~ 1s\deftined.) This 


luminosity is related to a change in magnitude Am by 


wwe 1070-4 4m (2.5) 


where Am is taken relative to the magnitude of the eclipsing 
binasy just berore the stare of the eclipse. “Ihe valuc of 
2 at the minimum of the eclipse is denoted by A. A quantity 


@ = a(k,;p) as also detined; representing the fractional 
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FIGURE 4, THE GEOMETRIC RELATION, 
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light loss. At any eclipse phase (see figure 2). This 


quantity may be determined directly from the observations 


by 
 Ll=-28 
Qa = (aoe (y 6 ) 
X = 1 - aL 


Zee lOcalvand Annular holupses 

We are now in a position to describe Russell's 
method for the determination of es on and i. The follow- 
ing derivation may also be found in a recent book by Kopal 
(HOT 9 apo. lO). The basic 1dea or the method 1s to write 
down the geometric relation (eqn. (2.4)) for three eclipse 
phases and to consider Shere. Bee 4, and ee as the unknowns. 


For the system of equations to have a unique solution, 


tee 2 
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SELie oe (1+ kp) i 


In this equation, k is the only unknown, Since p can 
in principle be determined from a(k,p). The phases ony and 
om) are Chosen so as to Correspond to Ga — 0-6 and a— 079 
beSpeCtIVvely. | HLOm Chis) pOlmte, on, che mecnod used storyde-— 
termine Cor ae and i depends on the type of eclipse. It 
should also be noted that the entire light curve is not 


Fequirederor the analysis... Only oneshali VoL san eclipses 


required. 
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Borsa total or annular eclipse (see figure 5), the 


determinant wan equatton (2.7), is written as 


sin’®, =i Bae sre (2.8) 


where 


A = ne ' B= A = sine, 


and D 3 
2(p,- P,) a k(p3 - py) 
W(k,p,a) = ae = soe ; (2.9) 


2 
2(p, DP.) + kK (pj) ~ Py) 
BE sin*e, is allowed to represent any eclipse phase, then 
A and B may be determined, and finally wW(k,p,a) for the 


given 9 Thus, one tabulates (k,p,a) for all eclipse 


3° 
phases. ~wtk,p,c) canealso be tabulatedyusing equation (2.9), 
SO a comparison between the observed and theoretical values 
of y(k,p,a) can be made, allowing k to be determined for 
each eclipse phase. Russell tabulated w(k,a) for both types 
of eclipse, but the most comprehensive tabulation was that 
Of Russell and Merrill (1952). A shorter and more useful 
set of tables (for x=0.5) was published by Irwin (1962). 

To summarize, one finds a value for k by determining wv (k,qa) 
from the light curve by equation (2.8), and by using these 
observed values of wW(k,a), along with the corresponding 
values of a determined by equation (2.6), to do inverse 
inberpolationesin a table of W(k,o), thereby producing a 
Gange Of values for .k., An average of the values of kK is 
taken, and this number, <k>, is then taken to be the 


VOrctece km in later Calculations. Thesinelinatiton i ena 
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where o, (k) and >. (k) are two auxiliary functions also 
tabulated by Russell and Merrill in the reference quoted 
above. The value of r. can now be found by using the defi- 
Majelony OL Ke 

The method of finding ter ee and i just described 
was modified by Russell and Merrill (1952) to use more 
points on the light curve during an eclipse. This is 
achieved by using three weighted means of siné ana, U(k,a)7, 


and. by defining, a new tunccion Rix,k) 


M, (sin*@] - M,[sin“é] M,(¥1-M,{v] 
ene ea Premera PaO =r aa aot 
M,[sin 6] -M,[sin Sak 82 ee 
M,(sin‘6] = Ay + Bihan | ; Jie ea ee 
2 


where the ess vl, t= 1,273, LS a weighted mean of Share FOr 
certain values of a, and eels J=2., 273, a5: the correspond— 


ing weighted mean of y(k,a). Only one table is required to 


ae eo” 
il >{sin 0], and M,[sin Ol. 


The values of toe Be and Lare Obtained as an the ecarlaecr 


fund k Given hc, ks), 01 sin-6], M 
version of the Russell method. The version just described, 
known as the 'Russell-Merrill' method, has the advantage of 
simplicity and greater computational speed, and will be used 


in subsequent examples. 
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In the case of a partial eclipse (see figure 5), a 
different approach is required. The problem is more diffi- 
cult to solve since observations of both eclipses are 
required for a unique solution, and the value of a at mid- 
eclipse (denoted by ao) is also unknown. At mid-eclipse, 


the geometric relation (2.4) becomes 


De 2 2 
= + = 
Cos i ru kp) ag p(k,a.) (22D) 


Since. = Cat thus pornt. Upon "subtracting Enis result 


from the geometric relation, one has 


io De ee. 
Sin ne, esc i)[2(p-p.) + k(p =e ; (213) 


This is the fundamental equation for the analysis of partial 
eclipses. Russell's approach (Russell, 1912b) was to first 
define n as the ratio of 1-% to 1-A at any eclipse phase, 
and to take the value of sin’é abvin= 005 as a “base pointe 
Russell then divided equation (2.13) by its counterpart at 


n= 025, Obtaining 


2 2 
eae BCye is ce ies) fel) 
Sin 8 (a) _ = ____9_______3- = x (ka 50) (Qa) 
San 6 (0).5) 2(p)-P,) + K(p) - Po) 
where © (0.5) “denotes “the value on Go when n= 055, “@(n)) the 


Value Of. 0 fOr any OLner nn, and Py thes value of (peat on —.0r 5. 
tn the vanalysiceor partial ~eclipses, scenes funce1ons splay a 
Similar role as do, the W-ftunceions in chewanalysis .of sectall 


and annular eclipses. “However, an the partial eclipse case, 
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Ehemwcolucron Ls graphical. The values of X (k,a ¢n) may be 
computed from the light curve by using the left-hand side 
Siecocuattvon 2.14). Since ke and a. are co be solved Eor 
first, only two values of y are needed. Suppose these 


values to be dencted by y=c. and eS Coe One may also com 


1 
Duce. y from the right=—hand side of equation (2.14), and 
therefore tabulate X(K,a 7m). The most complete tables of 
X (K,a07n) are those compiled by Russell and Merrill (1952), 
which can be used for both occultation and transit eclipses 
and any value of limb darkening x. After choosing a value 
Of X, One uses the y-tables to plot a, as a £UNCELON Of kK 
for the given values of y, namely or and Co: Therefore, 


the point at which these curves intersect should provide 


the required values of k and ao: Unforeunately, the soluq 


tion “Obtained ls indeterminate since Wb as noe Known whether 


the given eclipse is an occultation or a transit. Therefore, 
both eclipses must be used. The type of eclipse may now be 
determined quite easily by using the relationship (see 


Piewane (LI 62 Pao O07.) 


x°° (ka 7n = 0.8) > x" (k,a 7n= 0.8) (D415) 
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relation may be verified by consulting the appropriate 
tables "Ore. SLOr nN =)0.o.) sAnNOthen Pproblemarises Inpunemiact 
Ena tethne antensecereion. cm, thes two — Constante wcunves can be 
Guite shallow, resulting in an indeterminate solution once 


again. To remove this indeterminacy, another independent 
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relation must be introduced. 
If A denotes the value of minimum light for either 


eclipse, then 


A=1-aL L. + = 
5 ( s 1) 
where L is the relative luminosity of either star. If one 
writesthis out for both eclipses and solves for Os there 
results 
LA 
oc ; 
fo) =e Meo Omi 2 es for an occultation 
e) a 2 
~ Y 
(2.16) 
ce 2 : 
= - + _ ar 
aS z: dy was h)k x for a transit 
where do and AD represent i for occultation and transit 
: : n : ac eae. 
eclipses respectively and Y denotes the ratio a6 pes . 


Either of equations (2.16) is known as a "depth" equation, 
Since such equations relate the depth of an eclipse (1-)) 
to a, and k. Equations (2.16) are incorporated in the solu- 


tion method for partial eclipses by obtaining values of k 


and Los for successive values of Y. Tables of ¥(ao7k) exist 
(Irwin (1962) gives tables of re Ca pee es A see for 
this purpose. The set of values of k and a, so obtained are 


plotted on the same graph as are the equations for x = 
constant mentioned earlier. The curve so defined will 
usually make a steep intersection with the x = constant 
curves, thereby rendering the solution determinate. An 


example of such a graph is shown in figure 6. 
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FIGURE 6, THE GRAPHICAL SOLUTION FOR 


PARTIAL ECLIPSES (ADAPTED FROM IRWIN(L962,°6,607)), 
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The preceding paragraphs have described the form of 
the Russell-Merrill method that is most useful for the 
analysis of total and annular eclipses, but which for par- 
tial eclipses is not the best nor the most useful approach. 
In fact, Kopal (Kopal, 1979, pp. 113-114) has argued strongly 
against the use of the x-functions for determining the 
orbital elements. The essence of his criticism is that the 
position of any "fixed" or ‘base’ points; as used in the 
Russell-Merrill method, can be determined only to some 
Finite accuracy, and that this uncertainty would propagate 
through the entire light curve solution, leading to an un- 
7a and i. Furthermore, the 


solution as fitted only at the *fixed points*, not at all of 


certainty in the values of r 


the data points. One could imagine a "worst case" in which 
a rather large initial error would propagate and magnify 
through the solution, leading to wildly erroneousresults. 
It is situations such as these which have led other workers 


to use other versions of the Russell-Merrill method. 


2.4 A General Formulation 
The Russell-Merrill method may be restated in a form 
useful for any type of eclipse, and moreover, in a form that 


is amenable to use with electronic computers. This method, 
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due to Kopal (see Kopal (1979), pg. 115), takes the geometric 


relation and rewrites it in the form y=ax+b, which is 


linear. If one defines 


x = sin-e and y = (1+ kp) 
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then the geometric relation may be written in the form 


x = by, oe ieee - cotti 
or 
Stes cos*i 
r 
g 3) 


this latter form being suggested by Tabachnik (1973). 
The elements et eae and i may be found from the following 
equations, where a= sin*i/re and b = cos*i/re 

tan*i=2 , x Cpe? eee ne an 

g g 

If an initial value of k is used te determine p from a 
table of a(k,p), then the correct value of k will be the 
one that renders equation (2.17) a straight line. The 
straight line is fitted using the standard least-squares 


techniques. A good initial guess at k can be made in seve- 


ral ways. The simplest is to use the formula 
, ere" 
Sram 


where O° is che phase angle of first contact and 9” the 
phase angle of second contact (see figure 7 for definitions 
of @6' and 6"). Other methods are given in Appendix 1. The 
advantage 1m Uusingmequacion (2.17) I2es  1n che vtace thatsall 
available eclipse data are used, and no special points on 
the light curve are required. Moreover, one can use any 
Mike etable wfor either an OCcUulLtTalLion, Of ayuLans Ue, sand 
for any lamb Garkening te determine pl(k,0). “Therefore, this 


version of the Russell method is clearly the preferable one. 
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The=present author has written a program for the T1—-59* pro— 
Grammable calculator to use equation (2:17) in the analysis 
Ofeecclipsing binary light cunves. The values Of ance and 
Peake Used as input. “Thertl—59 calculator Vs pantreularly 
convenient since it has a built-in least-squares linear fit 
routine, which can be easily incorporated into a larger 
program. This program will be used in later sections when 
Pabuteular stars are considered. A listing of the program 
is presented in Appendix 1. A computer program, LINE, in- 
eCOrporating Tabachnik'’s method, is also listed in Appendix 1- 
Before discussing the application of the Russell 
method and the Russell-Merrill method to close eclipsing 
binary stars, it should be mentioned that several other 
versions of the Russell-Merrill method exist in the litera- 


ture: Most off these are due to Kopal, in partacular the 


iterative methods (based on equation (2.13)), which have 
proven to be very useful. These methods are conveniently 
summarized in the 1979 book by Kopal. A computer program 


incorporating an iterative method has been published by 
Jurkevich (1970). An important variation due to Kitamura 
(1965), which employs Fourier transforms of the light curve, 
will be considered in the next section. Some methods which 


are no longer in use are those of Scharbe (1925) and Fetlaar 


(1923). The latter method is summarized in a book by 
Psesevicehn GL97o)y watch also contains a descr pevoneoted 
method called the 'express method'. A recent revival of 


Kopal's iterative methods may be found in Look et al. (1978), 
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which also contains an interesting application of the depth 


equation. 


2.5 Non-Sphericity and Rectification 

Naturally, one cannot apply the Russell-Merrill 
method (or any one version of it) to all eclipsing binary 
Stars. Not all eclipsing binary stars have spherical compo- 
nent stars since there are inevitably tidal effects in any 
close system. Those eclipsing binaries with relatively 
short periods, less than about 3 days, will most certainly 
have some tidal distortion present, since the two stars 
involved will be quite close to one another (Kepler's har- 
monic law: P'«a ). There are often other associated 
effects. An obvious one is that one star will heat the 
other, the effect being a mutual one. When first discovered, 
this effect was called the "reflection effect", since it 
was believed at the time that light from one star was re- 
flecting off the surface of the other: Though inaccurate, 
the name stuck. As the theory of stellar atmospheres 
evolved beyond the well-known "gray" case, it was realized 
that the "reflection" effect was really a heating effect. 
The reflection effect has become one of the most difficult 
effects to understand, and hence model, since the magnitude 
of the effect depends not only on the closeness of the stars, 
but also on the state of their atmospheres. The problem as 
it currently stands is summarized by Sahade and Wood (1978). 
A comprehensive study of the reflection effect, typical of 


many done, is that done by Napier (1968). 
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Another effect present in eclipsing binaries is a 
direct consequence of the closeness of the component stars. 
thas 2s the presence of streams of matter between the stars: 
Matter streams can arise in two ways, the first being the 
presence of one star with a moving atmosphere (stellar wind). 
WolT-Rayet Stars and red Giant Stars can be involved in 
this type of mass exchange. A mass exchange can also arise 
if one star expands out to its Roche limit (see figure 8). 
Some of the matter from the expanding star is then drawn 
Onn by the oOtner Tsitar (through 1ts*grayitation)i,, wWLen ene 
resul Pp sbeingweither an sccretion disk or a “hot spot’, 
where the matter stream makes contact with the atmosphere 
of the attracting star. In the Russell model, the effects 
Of tidal interaction, reflection, and mass transfer are 
dealt with by the process of rectification. 

The dynamical and physical theory upon which the 
process of rectification rests will not be developed here. 
A comprehensive treatment may be found in an article by 
Martynow (1973) an the book edited by Tsesevich (1973) - 

The treatment to be followed here is that given by Proctor 
and Ginnel ll. (L972). ThesRussell model treats “he stars of 

a close eclipsing binary as prolate spheroids (see figure9), 
although the nresults obtained at the end o= Lhe Tectitica— 
tion process can be converted into results applying to a 
tebaxial ellipsoid. “The objeCcr or GecCelmication 15) co 
convert the light curve of an eclipsing binany consisting 


Of distorted stars into an equivalent “spherical” light 
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FIGURE Go. THE ROCHE SURFACE. 
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FIGURE J, A PROLATE SPHEROID, 
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GULVe~ ~ Thus, Lecliftication produces a rectified luminosity 
ae and a rectified phase oe, given unrectified values of 2 
anc 0 MOBS hIGOrously,. If sane Observer at point 0, having 
divreceLon cosines 2,m, nm is watching am ellipsoidal star sin 
a Close binary system, then the process of rectification is 
an attine transformation that) carries the observer at 0 to 
another point 0°, with direction cosines. @° 7m", n', ateche 
same distance from the sphere. Since the transformation 
carries an ellipsoid into a sphere (with a radius equal to 
the ellipsoid's semi-major axis), then the luminosity of the 
spherical star, as seen at 0', must be the same as that seen 
from the ellipsoid at 0. Since the light from the ellipsoid 
varies with phase, the light from the sphere must be modu- 
lated to produce the same light variation. The ‘affine 
transformation is illustrated in figure 10. It was shown by 
Russellvand Merrill (1952)—8 that the light variacion from. 2 
prolate spheroid with axes (a,b,b) is the same as that from 
avteGlaxial e€lliiosoid, having axes (a,b,c). If thesorbit has 
an inclination j in the prolate spheroid case and i in the 
ellipsoid case, then the semiaxes b and c of the triaxial 


ellipsoid are related by 


Ds, » 
ee eee (Oo) 


aac De 


Before proceeding further, it is necessary to define 
some parameters that will be used in the discussion that 


foblowse Nhe Ooblateness © 2S cCetined vas (a—b)/a, a and Db 
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triaxial ellipsoid. An approximate expression for ¢« is 
eee” ! oe 

2 | W being the eccentricity of the cross=section of the 
Stal in Che Orbital plane. Ti will also be useful to de— 
ae gee hae tae ee 

Pe 6a eee eS. Sars Lo 2 Sari he A non-Svpherical star 
will not have a uniform surface gravity, and consequently, 
tyOse Parts of the surface Gf the stax tarther from the 
Star's cent €r will appear cooler, while those parts ¢ 
PO Eie Genteer (iear the pole of rokaltion) wiil a 
hotter. A Quantity that describes this effect 1s the 


gravity darkening coefficient y, which is defined by 


| 


Marnyney (L973); 
G 
a Z (2520) 
ec, /xT., as 
AME CA “ee ) 
fe) 
where X is the wavelength of observation, T_ the surface 


Ceniperatbure, and’ c, (equal €o he/k) is a numerical “constant 


2 
whose value depends on the units of } and T_ (see Gray 
= fe) 
(2276), ba. 117). Therefore, the observed antensity at 
any point on the star's surface will be 
b= Mit = 2s cos (i= tyes) ARC 
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where H is the intensity at the centre of the o 
x and y are defined in equation (2.1), g is the surface 
Gravity at any Poaine on tne stan, and Ty is a reference 
Value Gf Gg, usGally taken tO Be the value at the pole of 
tne star. The lignk from Grther Star can De Sxpressed “as 
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R(6) = 190°) (1"= we sin74 cos-6) +c £4) (Ze) 
where 
GS eG ear’ 


Meas era eae Yh 


In this formula, 1(90°) is the light from one star at qua- 


drature phase (8 = 90° 


), G is an ‘albedo Eactor that deter— 
Mines the fraction of laght received from the companion 
star which is reradiated at the wavelength of observation 
(Rossetti and Mern i) (1952) 7 pg.46) prand: £(o)=is a ‘phase 
Function “that. characterizes the rerlection eftfect.. The 
procedure of rectification is one in which the eifects of 
reflection and ellipticity are removed by writing out the 
equation for %.(0) for both stars, finding the sum of these 
two equations, and then doing the appropriate subtraction 
and division to produce the value of £(8) for a system 


GOnSiSteang Cf Spherical starss, ine praccice, secre icacion 


is done by fitting a Fourier series of the form: 


Cc) A, +A, cos & +A, cos ora B, sin ¢ F B, Sin 28 G2vaad 3°) 


to the light curve an the eclipsing binary outside the 
eclipses (one may use cos“ and cues instead of cos 28 and 
Sim2e5 bY Making wSe Of ay Eregonometric 1denurty, Duc the 
coefficients will then take on different meanings). The 
series may be fitted either by the least-squares method, or 
by a graphical method developed by Russell and Merrill. An 
example of the latter may be found in Appendix 1. The-rec= 
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and the rectified phase by 


Satie , 
sin 9 oa oo. 
Beet (l-zco : ts) )% 
: obs 
and 
*% 
Lis SZ 
Cos ¢ = cos 6 A eo eee Gy) 
rect eer | obs 
obs 


where both equations are required for proper quadrant de- 
finition, (in taking an inverse tangent, an. electronic 
COMpucer USes the range =1/2 <9 "< 7/2, rather than O< 6s 2, 


which is the range of os ). These formulae apply to all 


bs 
Observations, DOth in tand out of eclipse. The factor Zz may 
be obtained empirically by using the approximate relation 
fe | 2251 - It Should also be noted that the presence of the 
Sine terms in the Fourier series for 2(6) is not justifiable 
physically; their only purpose is to take care of any extra 
VoomplicacLions. Chat Might arise.» This then 1s the precess 
of rectification as developed by Russell and Merrill. The 
end sprOduCcte Seas LigGhtwecunve= chat sis elat outside the 
eclipses, with the eclipses being those appropriate to 
spherical stars. 

The process Of FeECtLElcation 2S Open co criticism 
On Several grounds. The first, and most Obvious, isthe 
presence of the sine terms in the Fourler series for &(90). 
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SUGN  JUStification, published or unpublished. ~Another cri- 
ELGUSM;p catsed ‘by Kooal (1979, pds 192) > ls that onesls using 
@ePOuUbLer Series outside its range ot validity, since 4 
series, which has been fitted to the out-of-eclipse obser- 
vations is being applied to all observations, both in- and 
SuEsOL—-eclipse. PRectriicatzongwikl workeator systems! in 
which distortion effects and the reflection effects are 
Minimal. The example in Appendix 1 is of this variety. In 
cases such as these, the B terms are quite small in com- 
parison to the A terms. However, one really cannot apply 
rectification to highly distorted systems (very close 
binaries, e.g. W Ursae Majoris). In systems such as these, 
the shapes of the stars depart greatly from an ellipsoidal 
Loum, and actually approach) avy Roche-=surrace form... The 
theory upon which rectification rests is clearly not de- 
Signed with such systems in mind. Consequently, rectifi- 
cation is no longer used, and more physically acceptable 


procedures have replaced it. 


26 J DLEterential, Corrections 

If the geometric elements i rly, i, Ly rl, , Xp 
and x, are well-determined (in the sense that the solution 
for these elements is determinate), one may improve the 
values of these elements by the 'differential corrections’ 
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then 
AZ(§8) = Au —- aAL - LAa 
N 5 
=Au-aAL-L ) ax. , (Ooo) 
noe x, J 
ue 5 
where x. 1 i ] ; 
: s one of the elements Yur for, Xr ene X5 


Equation (2.28) may now be regarded as an equation of con- 
dition, Soe thatelfsthiss equations written out for teach 
(9,2) pair, one may solve the system of equations for Pa aN 
Au, and AL by the least-squares method. The value of A is 
found by subtracting the calculated value of 2 from the 
observed value (i.e., an 'O-c'). The various partial de- 
rivatives appearing in equation (2.28) have different forms 
according to the eclipse type. The form of a must also be 
chosen according to the eclipse type. The paper by Irwin 
(1947) describes the procedure of differential corrections 
in great detail, and tables of the various derivatives are 
provided in an appendix to the paper. The present author 
has written a number of computer programs for performing 
the differential corrections procedure using the values of 
the derivatives from Irwin's tables or values generated by 
the equations for the derivatives. Some of these programs 
may be found in) Appendix J. Dt should be noted thacvone 
Cannot apply differential corrections to vevery eclapsing 
binary star, since, as mentioned earlier, a well-determined 
set of elements is required, as well as a large number of 
observations to make the least-squares method truly appli- 


cable. Least-squares differential corrections should not 
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be regarded as a 'black box' that always generates improved 
values of the elements, therefore it should be applied with 
some discretion. As Irwin mentions in the paper quoted 


earlier, least-squares is no substitute for good sense! 


2) CONCLUSION 

The discussion of the Russell model and the method 
of light curve analysis associated with it is now complete. 
This model of an eclipsing binary star is best applied to 
systems having spherical components, since any application 
to systems having distorted stars will inevitably lead to 
the use vot rectification, the validity of which is in some 
doubt. The Russell-Merrill method is still used to provide 
a preliminary set of elements to be used in more advanced 
methods of light curve analysis. In short, the Russell-— 
Merraliemethod 2S not tehe most, ©£ruLELUL One, since ae ws 
possible to derive much more information from a light curve. 
It would also be of some advantage to have a method of 


analysis tailored for use on an electronic computer. 
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CHAPTER 3 
THE METHOD OF KITAMURA 


S20. cintroduction 

With the development of the electronic computer in 
the lace 1950s) many workers in the iield of close binary 
stars began looking for methods of analyzing eclipsing 
binary Light curves’ that might be suztable ton electronic 
computation. The programs that were developed all used the 
Russell model and the Russell method for computing the 
orbital elements. In general, the method of analysis used 
was one of Kopal's iterative methods. A good example of 
such a computer program can be found in the work of 
Jupkevicn (£970)=— Thais’ preliminary analysis was asually 
followed by the application of a differential corrections 
program to obtain improved values of the elements. However, 
the problems of determinacy were still present (especially 
for partial eclipses), as well as the ever-present problem 
of allowing for the various ‘proximity effects' present in 
close eclipsing binary stars. It was suggested by Kopal in 
1959, that the problem of the determination of the geome- 
tric elements of an eclipsing binary might be more easily 
solved if one were to somehow make use of the Fourier 
transform of the light curve. Such a method was formulated 
Dyk camuras(1965). selce stills usedmas an DasiomonenRusseial 
model, of Jam eclipsing binary, with Lect iication being 


required fon close eclipsing binaries. | Kitamura s method 
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solved >in pare the problem of “proximity effects", since 


the Fourier transform operation (like any form of inte-— 


Gratton) 1S a smoothing operation. | The required Fourier 
transforms being easily calculated by an electronic 
computer. The process of determining the elements was then 


one Of matching the Foulrier transforms of the observed light 
curve with those computed from the Russell model. This is 
the essence of Kitamura's method. We now examine the 


method in detail. 


3.2 Incomplete Fourier Transforms 

As previously mentioned, Kitamura's’ method is based 
on the Russell model. The observed light from an eclipsing 
binary system at any phase 9 is, according to the Russell 
model, 


VCO) = i= ola 


where £, a, and L are the quantities defined in chapter one. 


Following Kitamura (1965, og. o0)l, the Fourier etranstonrm of 


Separating the real and imaginary parts, we arrive at the 


‘incomplete' transforms of 2 (8) 
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BOr Che purpose Of analyzing the light curve, the Lert hand 
Sides ob “equations (3.1) are Used) white the sight Rand sides 
ane Used ton the model computations. §The specification of 
the limits of integration remains to be decided. Since 
£(8@) is symmetric about either minimum, the lower limit of 
integration is obviously @=0 (assume that @=0 specifies 
the minimum of 2£(8@)). To determine the upper limit, one 

8 


may choose any phase o> 6 being the phase angle of 


Ad Ik 
last- contact (see figure ll). The eftect of choosing «@ 
arbitrarily will be accounted for later. It should be 


noted that in general, there will not be a data point at 


exactly 0 = 0, so any data point having a phase <« close’ te 


8=0 will suffice. Equations (3.1) become: 
b b 

2a = | 2 (a) Sin, mesd9 and c= | Vai OVeCOsSenecds . (3:72) 
] 
€ E 


Without going into excessive detail, Kitamura found 16 con— 
venient to define three parameters simply related to Sy and 


Co to connect the Fourier transforms to the geometric ele- 


ments. These parameters are 
lh. S lb S I= © 
Loe: bee Lo 
F — 1 F = and E = = (hoe) 
al Lic, 2 Ljc, 1-A 
where 
lec. =o = rh =~ Cc ; (Gas) °=.0"=) 
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FIGURE 11, THE ANGLES ©,9),9, 
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Once again, X is the value of 2(@) at mid-eclipse, and K is 
esnoLtmalization factor for the light curve, which) can be 


computed by 
taney 


or and >. being two phase angles outside the eclipse. The 
birse two) terms, 1m each o£ equations (3.4) account for the 
choice of integration limits, ¢ and e«. The parameters Poe 
Foo and E may also be calculated with the right hand sides 
Of sequations (3 j= Kitamura <¢1967)| pextormed jsuch <a com-— 
putation, and the result was a large set of tables of Fie 


Foy evqrel J Smena “OP 6G Ser <a sane AO so eh BS Oe TOmaGCcOun tutor 


pe 
limb darkening, Kitamura computed Foy Fo, and E for zero and 
total lamb darkening. jHe then constructed) a table of “re 

lated delta-functions', which allow interpolation for inter- 


mediate values of limb darkening. The ‘related delta- 


functions' are defined by the following formulae: 
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where a prime denotes an occultation eclipse and a double 
prime a transit eclipse. For any darkening x, the charac- 


teristic functions become 
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£(0) being the light loss at mid-eclipse. If the left hand 
Sides of these equations are taken to be the values found 
via equations (3.3), then if preliminary values of Pi rts 
and i are known from an initial solution, the values of the 
A-functions may be found from the tables and (assuming a 
limb darkening x) used in equations (3.6) to solve for a 
Sal and Eo: Using these new values of the characteristic 


functions, one may find new values of tis roy atcieias 


geo Practical, Approach 
In practice, the incomplete Fourier transforms oa 
and Co are conputed from the daght curve by using the tra— 


pezoidal rule: 
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Le obtain Lmproved) values vor So and Cae end corrections to 


the trapezoidal rule may be used: 
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a) 2 ini | 
AC = 0 : ACL =D We )iecinn ome AC, =<% ho )rsin 2o |; 
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AS MN A oe ye = Zh*[2(0) - 2($) cos 26) 


ee 
Teerot 
With these values of S and C , the values of L.s and L.c 
n n Lee aie 
in equations (3.4) may be found, and hence the values of For 
Foy and E. The preceding analysis is performed for each 
observed eclipse, provided that the eclipses are deep 


enough to allow Sy and co to be computed. If the eclipse 


type is in doubt, the following 'rule of thumb' may be used: 


This rule need only be applied in the case of partial 
eclipses, when the type of eclipse cannot be determined 
divecely £rom the light curve. “One can now consult 
Kitamura's tables of Fie Foo E and FL/E, LOuetind Ene taper o- 


priate elements Pir tos and i. In general, one will end up 


2 
With ceveral sets of elemenus tor which the: Values of the 
Characteristic functions match. To determine which set of 


elements is the ‘best £1t", one can examine the values of 


T-1, where T is given by: 
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where the values of as and ie can be computed from the 


EtCSot Ot cauatlLonse (3.4). and Eee and a Ue can be found 
from Kitamura's tables (ceo SS Pere ola cnet = OGG, in. ane 
fo) a O O a 


Q are listed in the tables adjacent to each set of elements 
Tirta, and i). The set of elements for which T-1l is a 
minimum (= 0) is the 'best' set of elements. If limb dar- 
kening is to be taken into account, the procedure described 
earlier, involving the related A-functions, may be used. 
The present author has written a computer program 'LCFT2', 
which computes a er Fie Foy 


eclipse light ‘eurvye> "fhis program, along with 2 saniple= run, 


FL/F, and E for’ a-given 


may be found in Appendix 2. 

As with the Russell method, a differential corrector 
pregram may be used. to improve the values of the bese set 
of elements. fn face, aM dlrrerential corrector sheulaepe 
used since the values of the elements obtained are read 
directly from the table and therefore represent only an 
approximate set of elements. The programs for differential 
corrections discussed in chapter one may be used for this 


purpose. 


324 Conclusions 
Kitamura's method has many advantages over Russell's 


method, since it uses all the data available instead of only 


43 


od most ino ipaes 6d aes | 

bad sth ae ho)" Ge” 19) a iolenen 

hen .5 faxed = OF ror a 5 es > a ext @ 
as 

dinteseia ta s9@ doee, od sachin aol ods mi Bowath-o 

it dotdwesot asso se soa oT ~(h Bae Pug 


& a2. 
-xeb dmkt tL adnsinals 14 t98 ae ared et (0 =} 


swuiiscveh siubpoezq Sclu {MEO DOE cyst Hate od Of Bi | niaok 


ed wan  Gqavioneae eats ent pnéviovad aehtas 


e 
MOE tt 


‘UTA’ morPeTa astranes. G pedelaw 26d xiitos sepeete 


4 = wi. 


qavle 4 moa a. Sav ay oe q vy) Ped " Pd aed OgMaS 


(met Simgs 6 ftw. trois _yiivapes audit .aedtiy, dnobi — 
.< xibieogs nb bane? sd-4 


(arene tests 8 © \ borden Llegene eat Rete ef 
. : coe ae 


tosbes20- 
453 saad off Se barley ‘Sst avenge 64! ieee scl yan 
Sd bivode Ao): i ~~ — n naar. 5 «ines An nsoaagte | 


GLes SSe nopladds ant ai] es oes aa heviiev ad sone 
a 
c tb + 


7 


5 L, 
iv 


a6 vino digeamqes Sxetingst! his eles eid meek * 
Pala ial tf 203 Kine ®So2q st? .2fnemeie 29 yor oter 


sido vot Seed on (an eto verqadt me hoeauoa:l 3 


44 


three or more fixed points. Furthermore, as mentioned 
earlier, the Fourier transformveperation "acts vas “a ssmooth— 
ing operation to reduce the effect of random errors in the 
observations. This aspect of Kitamura's method also has 
the advantage of speed and computational simplicity. Onee 
the incomplete Fourier transforms have been found, it is 
just a matter of consulting the appropriate tables to get 
the values of the elements. The bulk of the computation is 
done by the computer. 

However, there are distinct disadvantages in 
Kitamura's method. The first is that it is based on the 
Russell model of an eclipsing binary, so FECtliacaclon as 
required for those cases in which the stars involved are 
Significantly distorted, and where other effects (reflec- 
tion, mass transfer) are present. Kitamura's method would 
be much more useful if it were based on a better (more 
realistic) model of an eclipsing binary. A second disad- 
vantage is the eed cS consult a rather large set of tables 
in chery dast Step. im the analysis cf a Light curve. This 
procedure is not objectionable when results from both 
Minima are available, or when a single minimum provides 
determinate results. A search through the tables in cases 
such as this takes only five or ten minutes at the most. 
However, in an indeterminate case, the researcher is faced 
with the prospect of wandering through 241 large-format 
(about 30 cm =< 50 cm) pages to look for one set of elements. 


One can easily spend an afternoon doing this. All objections 
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aside, Kitamura's method provides an excellent way in which 
to compute the geometric elements of a well-detached (no 
photometric complications) eclipsing binary. From a his- 
torical point of view, Kitamura's method provides a start- 
Ing pOmnt Lon an alternative formpor Pught curve analysis, 


which culminated in the recent work of Kopal (1979). 


45 


ahi emote S HO 


7 
aid sme ae Let iios yesecy 
sabivoxe borden s' swiss LP «week ‘40 sab09 


sia be 
= - we 
ehey Lams sytyn, dnp iy yo sriot Svitenias (e “ws 20% miLog & 


~* 428 = ; 
Jzace" e - 


E 


- 
(eee) Teger Gu grow sasces ed! ai tesantulus fae J 


CHAPTER 4 


ANALYSIS OF LIGHT CURVES IN THE FREQUENCY DOMAIN 


Avie unerOoductilon 

To make further progress in the analysis of the 
light curves of eclipsing binary stars, one must dispose of 
the Russell model entirely, and adopt a more realistic, and 
consequently a more physically complex model of a binary 
Stan. Brierly, such a model would take into account details 
such as non-sphericity (e.g., using the Roche model), stellar 
atmospheres (for the purposes of deducing limb darkening, 
gravity darkening), and any other effects that might be 
present. Consequently, the mathematical machinery required 
grows in complexity. 

A generalized approach, such as the one just des- 
cribed, has been developed by Kopal. A convenient summary 
Of this work cam be found in has recent, book (1979). Briet= 
ly, Kopal considers eclipsing binaries with both spherical 
and non-spherical component stars from both a dynamical and 
a geometric (i.e., eclipse geometry) point of view. One 
can, therefore, look at an eclipsing binary in any degree 
of complexity, from a simple preliminary analysis to one 
takinge into account thevunigque, teatures on a particular 
binary star. Kopal achieves this flexibility by using the 
PourLter transrorm of the lightvcurve, but gin alway chat as 
MUCH GL Eferent, Lrom the approach used by Kitamura. | The 


details of Kopal's use of the Fourier transform of the 
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light curve will be expounded later. The motivation behind 
this approach 1S cCwoLola thes ourier enalysis, 1m certain 
cases, renders the solution for the geometric elements 


(ry, tr, and i) algebraic, in the sense that the end result 


ri oss, LOr example. Furthermore, 


"proximity' effects present in close binaries can be separ- 


1s an equation stating 


ated from light variations due to eclipses by the nature of 
their frequency spectra. Eclipse variations produce a con- 
tinuous frequency spectrum, while proximity effects produce 
only discrete frequencies. In short, Kopal's approach pro- 
vides both a realistic description of an eclipsing binary, 

and a relatively simple method for determining the geometric 


elements, given the Fourier transform of the light curve. 


4.2 The Equations of Kopal's Method 


TO begin, oun exploration of Kopal"s method, let us 


consider the Fourier transform of a light curve. The pro- 
cedure and notation follow that of Kopal (1979). Consider 
fiest the buch, curve: drawn an Petact es (il glen peer) 
coordinates. The motivation for doing this will become 


clear later. Figure 12 shows the light curve plotted in the 
v= eae coordinates. The case of spherical stars will be 
considered first. There is a straightforward generalization 
for non-spherical stars. Let the area bounded by the light 


curve and the lines singe = 0, & = 1 be denoted by Aan 


This area, known as the m-th moment of the light curve, is 


evidently given by 
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FIGURE 12, THE DEFINITION OF Aon, 
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Mon | (1-2) d(sin 


Se eR (pe Ie ee (4.1) 


where oy is the phase angle of last contact. The upper 


en: te ay was chosen because 1-2 equals zero for 6 > Os 


Giving Zero COnteibution to the integral. By the very 


nature of the integral defining A , one might suspect 


2m 


that Aan To Telaced te the FoOuLrLter tleansLtornm of L— 7 (6)> 


indeed, Le 1S Mot Gitricult=to show this. To preceed, the 


EGLIOWwinNG Lesult 21s required (Oberhettinger ,19/3; p2 31)% 


5 
M ek) tacos 24.0 
shqrilg Ra a (4.2) 


qm 0 (m+)! Gn=7) ! 
where 
i j = 0 
es A 
J 2 ano 


The differential of this series is, after a little algebra: 


(lye rein 46 


mtj)t(m-j)! eee 


Mulgiolying by l=.) integrating Leom 0") co Oo. and inter— 


changing the order of summation and integration on the 


right-hand side, we arrive at 


8 
Ay = : Tey Baa 
ae (1-2) d (sin ) 
0 
ay Coie a 
= 


To relate this result to the Fourier transform of (1-2), 


Gensidger anearbatrary function £(0)e(whieh comdsbesl-0(0)), 
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whieh 1S an even functron (isce., £(-0) =—-£(e)). The 


Fourier transform of £(6') will be 


Cc 
riv) = | £(@) e7272% ae, iS) 


where v is the spectral frequency, and tc are the limits 
Of Integration, which can be Gntinite., In the situation 
at hand, the limits of integration will be the phases of 
first andsltast contact sance £(6) = t-((e), which will be 
zero outside the eclipses (/6§|}>c). Let us now split up £(96) 


in the following way: 


a) 

fap) 

I 
|r 


[£(8) + £(-8)] + 5[£(8) - £(-8)] . 


The first term of this equation is an even function, while 
the second teym is odd. Knowing this, 1 as not dirtiiculte 


to show that 


(e Sc 
| 5 [£(8) + £(-6) ] Be ee ays = | [£(6)+£(-8)]cos 2nrvé6 ae 
a 0 
=F, (y) 
and (4 56) 
(e (eS 
| 5 [£(8) Deo eo | [f£(9) -£(-6)]sin 2rve dé 
S 1 ky 


The Fourier transform F(v) can then be decomposed into its 


real and imaginary parts as Fi (v) ~1F,(v). Now assume that 


£(=6) = —-f/(@), so that Bit) is zero and (letting c= 85 = 
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phase of last contact) 


a 
F,(v) as) | £ (QO )esine2inverce F (4257) 
0 
BOL (ey) = =e) and \) = 4/1 ene has for Foy): 
eae 72 eae 
5 EF. (=) = | (T= \ssimecperao =. (4.8) 
0 


This integral is identical to the integral in the series 


expansion for A derived earlier, and is also the expres- 


2m 
Sion £or the Fourier coertrticient ae in the Fourier sine 
series 
1-% = ay b,, sin 236 , (4.9) 


This series is a special case of the more generalized 


Fouriwer Series 
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ny! A 
b = 2 F (J) = = (l= 0) sanv276 de. (4.11) 
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0 
Banmadly, 
m aed 
eA nm (=1 iI 
1s ae nT ay mene ne (4,12) 


This result shows the relationship that exists between the 


pat 


oy enue 


: 


2G). Got see SAS ESS Dae sb 


(a) : aay died : 


a 
. - 


epiney alo, oh Lerponns ef3- om per al saga 6 


“gigs sft cela st ishe heagee fee Er Sep eM 108 4 
4 dhakgit2eae raiiont eft = 


ania selzagl SAT. mt: » 


4, 


| 7 
LB. hy Sponie. dt { = fri 
2 oe ee 


066.8) 
vebrstes ef (11,4 ni nolgeapezad to Simhl segge atts & 
yi seve cet> 225 | od -sneinstawe Fetes’ ser. 48% 


as 
an 
bate 


,. ab ort nieierad) a be cbs ate ae 


(ot 


pes a ae knit 


Ann and. thet Fourier: transtorms.o£f the light? curve..< The 
moments of the light curve, the A, Sr are essentially 
‘weighted means' of the Fourier coefficients at Formche 
CaseseOrcvinterest ansdlightecurvel analysis, m=152,3. The 
Aj, § are then 

A, = 7, 

eer Meg wore pe eee are Ue Ts (4.13) 
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A point of particular importance here is that the moments 
Ayn are related to the coefficients oa and these corres- 
pond to discrete frequencies in the Fourier sine series 
(4.9). The rmmdesxm on the Ao, 5 will then correspond to 
the frequencies used in the spectral analysis of the light 
curve. 

Another feature worthy of note is that equation (4.1) 
is valid only for spherical ‘stars. Furthermore, it applies 
for any value of limb darkening and any eclipse type. The 
computation of the A, functions will be considered later. 

The results just obtained may be generalized to the 
Case in Which ‘the Stars constituting the eclipsing, binary 
are non-spherical. From the point of view of data analysis, 
this step is quite simple. One need only replace the 
rnteqrandedrequacten (4a) (thang 1s, sthe fLuncceron—2(¢)) 
by @(n/72) — v(6) , where “(7/7 2)) is the value of (9) at qua— 


dvature. In the case of spherical stars,) 2772) = (0. )—1, 
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the upper limit of integration being 6 In the case at 


L° 
hand; however, the upper Jimit of integration is extended 
tOeU = I/2 since. (m/2) i= £9) 1s;8inegencral, non-zero for 

OM s50s= 1/220 hur thermore, esince at 1seofteny ditiicult £0 
separate the eclipse phases from the non-eclipse phases in 
non-spherical eclipsing binaries, an extension of the upper 


limit of integration to 9=7/2 will adequately cover any 


eclipse phases (between 9=0 and 9=7/2) that might be pre- 


sent. The integral defining Aan becomes: 
sh 
aT om 
ee | [2 (5) P(ejeld (San =o) -. (4234.4) 
0 


For future reference, this particular integral will be 


denoted by A to distinguish it from its spherical geome- 


2m 
try counterpart,” equation® (4).)* 

To implement Kopal's method, one must be able to 
determine Aon (ays Boe (whichever is appropriate) from the 


observed light curve of an eclipsing binary. There are 


OF Ale ;o Given 


several methods available for computing A om 


2m 
We) 2. 

The most Straightiorward approach ws to numerically 
evaluate equations (4.1) and (4.14) by the trapezoidal rule. 


If one samples the light curve at equally-spaced points, 


the standard trapezoidal rule applies: 
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where N is the number of observations, (8. , I-k,) are the 
Nedatas points, and = 9,/N is the sampling interval for 
the integral. Alternatively, the date may be used directly 


in the following integration rule (Niarchos, 1981) 
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Similar formulae hold for the evaluation of va The pro- 
cess of numerical integration just described is very useful 
and is adequate for most purposes. However, there are some 
drawbacks. First; One as Unable tO pertonm any sort of 
error analysis if the moments of the light curve are formed 
by numerical integration. The various aspects of error 
analysitc winkiebe exploned later. Moreover, om the case of 
the De a lack of points near 6 = T/2 will seriously 
affect the value of the integral. This effect has been 
noticed by Koul and Abhyankar (1982) and also by the present 
author in earlier work (1983). This problem arises from 
the fact that the value of 2(8) closest to @ = 7/2 gives 


small but significant contributions to A Since 2(6) is 


ue 
increasing in the range 0 < 8 < T/2, one would expect a 


large contribution to A a from values of £(6) near @ = 0, 


2 
but steadily smaller contributivons as @ approaches 7/2. 
This problem as Tllustrated in Licune 13. his ettece sas 
Particularly prevalent when’ am observer has concentrated on 
acquiring a large number of eclipse observations, but has 


paid relatively little attention to the Pight curve outside 


the eclipses. 
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4535  ThesUse: of Fourier Serres 

A much more suitable method for determining the 
MOmMentsS Of Ehewillght curve is to fit a Fourier cosine 
Sepiles to 1=0(0) or 2imn/72)—2(e) = In “the spherical star 


Gase, One ‘can fit: 1l-< directly with a Fourier series of 


the form 
N 
eg Wee cos (BT) (Ary) 
Ae 8) n S) 
n=L1 de 
where 6, = phase of last contact and N is the maximum number 


i 


of terms in the series to be used. The cosine series is 
used because 1-2(8) is an even function of 9, hence all of 
the sine coefficients will be zero. To find the various 
Pour ver coefficients an (y= 07 ee. NN) En eGuide VOnes aay 30, 


one must first estimate 6 and then form equation (4-17) 


i! ’ 


for each observation in the range 0<6<6 The resulting 


L° 
system of equations may be solved by the least-squares 
method. This process for determining the Fourier coeffi- 
cients is described in some detail by Kopal (1982a). Since 
the least-squares solution can also produce estimates of 
the uncertainties of the Fourier coefficients, an error 
analysis may be done. Before venturing into this, the 
relationship between the Fourier coefficients and the mo- 


ments of the light curve A should be explored. Following 


2m 
the development given by Kopal (1979, pg. 241), let us begin 
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AS yWwle equation (454) "or this chapter, ditterentiate both 
Sides OF this equation, multiply by l—=%, and integrate 


between €@=0 and §6=6.: 
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Aon = | (1-2)a(sin*™ 9) 
(4.18) 
= is) } eee aera (1-2) a[sin(2j5+1)6] 
2 10 P (m+j+5) 1 (m-j+5) 


Kopal has evaluated the integral on the right-hand side of 
equation (4.18) in terms of the Fourier coefficients a, 
The procedure used is not unlike that used to derive 
equation (4.12) earlier. Without going into the algebraic 


details, the following result is obtained (Kopal (1979, 
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The number of terms to be used in the second sum in equation 
(4.19) will depend upon how many coefficients can be deter- 
Mined sioqnaticangly from the data, Typically, one can 
determine only five or six coefficients from the least- 
squares solution of equation (4717) ] To faciulitate the 
determination of the Aan from the Fourver coefficients awe 
Kopal (1982a, pp. 131-132) summed the first series in equa- 
ELOMN(4...9) explicitiy, leading Lo much mores tractable 
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The present author has adapted a program published by 
Jurkeviceh (196 1), whieh performs the Fourier analysis GE 
the light curve and the computation of the moments. 
Examples of this analysis will appear in later chapters, 
when individual stars are considered. The process of 
least-squares Pourrer analysis 1s quite surticient for the 
analyeis Of most eclipsing binary Wsitars; (but Le mustcbe 
Genera ized aie tow to be appivedite verygcetase Gclapsing 
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W Ursae Majoris-type systems). This generalized form of 
Fourier analysis has been developed by Kopal (1982b). The 
mathematical development of Kopal's generalized Fourier 
analysis 1S quite complicated, and only the results will be 
stated here. The method also employs the dynamical theory 
Obeclose binary Stars, which will be considered Tater in 
this chapter. One begins by writing 1-2 as a Fourier cosine 


series valid in the range -7r/2 < 6 < 7/2: 


il M 
1-%== } e€ a cos2me , (Am) 
oe m m 
m=0 
where Em is the factor defined: in equation (4.520), and M is 
the number of coetficients to be used. Once again, the 


least-squares method is used to determine the Fourier co- 


efficients an The next step in the analysis is to deter- 


mine the ‘modulated moments' Bue from 
(Uy ee ee en) 
B =e ae ae COSTES SO oes CCN es 
n m m 
m=0 
The p*')_ functions have been tabulated by Kopal (1982b, 
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pp. 447-448). The modulated moments ee 


may now be related 


to the moments of the light curve through 


where the 04°") _ functions have been tabulated by Kopal 
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equations in N unknowns). Sance only Ags Ayr and AG are 
required, only three equations are required. However, AG 
is still to be determined. AG is now given by 
J 
ease an Cz or (4.24) 
j=l 
where 
= Ao 2 Clie apt Mtele Mette! el echt sy rove) tanec 
2° 0 a 2 M 


The coefficients Ss are part of “the Ligne curve modulation 
process, and will be considered in more detail later. A 


further equation incorporating the e5°& “Ss 


(4.25) 


Equation (4.25) may be used to solve for the ra hy matrix 


(r) 


a may be found from equatien (4.22), 


inversion), since B 


and since the second summation in equation (4.25) contains 


(A, 0) 


now be evaluated via equation (2.24). This completes the 


which cans berydetermineds tor, y= 1,.2),3.. oe can 
determination of the moments of the light curve for the non- 
spherical case. The procedure just described is a 'filter' 
£Or removing proximity etfects, namely ellipticity (non— 
Sphericity) and rerlection.” Tf tidal distortion, 2s i1mpor— 
tant, further results from the dynamical theory of close 
binary stars may be used. As a final note, one may deter- 
Mine the error in the resulting values of the Aan found 


via the process of generalized Fourier analysis. This is 


done with the aid of the following equation: 
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j=n J J 
where, once again, the ES a functions have been tabulated 
by Kopal) (1982b, pp. 445-446). Equation (4.26) gives the 


error Wncurred by Using n= 374 im equations (4.22) and (4223). 


(r) 


The values of the En "s will be quite small in most cases. 
4.4 The Use of Numerical Integration 

The moments of the light curves of non-spherical 
eclipsing binary stars may also be found by the technique 
of numerical integration (Kopal, 1979, p.195) discussed 
GCariierssRecalil that (equation (4.14) ), 


M2 
= Ee 


os T) - 2(6)]a(sin*™ 6) 


was the expression for the moment of the light curve in the 
non-spherical case. As in the Fourier analysis technique, 
one may modulate the light curve and obtain the moments Aon 
of the Light curve’ for “an *equivatent “spherical” system. 


To do this, the following equation must be used: 
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where B(m, 5 j+1) is the beta-function (equal to Pr (m)P (5 §+1)/ 
WG tee eB are the dynamical coefficients mentioned 
earlier. The last term in (4.27), the Bon! is known as a 
"photometric perturbation", and describes the departure of 


the star's shape from the spherical form. These quantities 
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are usually defined by rather complex expressions. The 
general form of the Bon ECEM alae (nODals, oL97 9, to. 195): 
2 
A+1 Ws 
See La ee, ee fie ee aaa 6) (4.28) 
2m it Hes * a ® 
< 0 


where Ly is the relative luminosity of the star being 


eclipsed, A is the degree of limb darkening, and the co- 


efficients es are defined by: 
MODs Hee ieee se 
; LU, 
1 - 
Ct +2 
and (A329) 
Pa) ae “h 
i re 
oe tae 
al eral 
the u,'s the cCoeckiicients of limb darkening (u, = linear 
darkening, u, = quadratic darkening, etc.). The term in 


2 
square brackets in equation (4.28) describes the effects of 


rOotvatltonal and tidal <istortion. To determine the Ajn 


for m= ,2,3), sone computes the Aes by Mumerical integra— 


2m 
tion, and the coeiiicirents os by numerical integration via 
a 
(a,n) 


= ae 
ahi ) 2) 2(8)]p, 
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"modulating polynomial', which describes a star's surface 
distortion to first order. Kopal presents these polynomials 
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(1979, pp. 198-203). Alternatively, one may evaluate the 


Sumernegquation .(4.27)) directly by.using 


a 
n 
ee 
m ) B(m, =jt+l)c. = [2 (5) Oe seeds ; (ey) 
= 2 a 2 m 
j=l es 
where the aes are polynomials similar to the oles 


polynomials. At this point, one would evaluate the Bo 


term using an expression appropriate for the eclipse type. 


Kopal gives B for m=1,2,3 for the case of a total eclipse 


2m 


(1979, p. 211). However, complete expressions for all types 
of eclipses may be found in the pages by Livanion (1977) 


and Rovithis-Livanion (1979). With the Boy term evaluated, 


Ene wpEeocess or finding A is one of simple algebra. Once 


2m 


again, it should be emphasized that this particular method 


@f determining the A, "s is useful only if the light curve 


2m 


has been well-observed, particularly around 9 = 7/2. 


{7 COMPU nGa the AO with the Kalman Filter 


An approach which is entirely different from either 
the Fourier analysis or numerical integration techniques 


for finding the A Ge is the Kalman filter algorithm. This 


eZ 
approach was first formulated by wWunkevich) @1976) in 
essence, the Kalman filter provides a means of determining 
the values of the Aon recursively, while producing a result 
that is optimal in the least-squares sense. It combines 


the features Of numerical antegration, and least-squares 
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the Aon 8 The theory of Kalman filtering has been inves- 
Cigetecdivin, detail, sand a particulerdy good <discussion uiis 
Provided im the book by Bryson gand Ho (1969). 

In order to use the Kalman filter, one must be able 
to write the equations of a problem as ordinary differen- 
tial equations. These differential equations are then 
treated as ‘stochastic differential equations', or equations 
in variables which are statistical in nature (random vari- 
ables having a variance o* and a mean equal to zero). These 
differential equations are then integrated between the 
observations, resulting in a set of difference equations 
for the variables involved. In the problem at hand, the 
variables are Am (hereafter referred to as Wr X= lek, 
and the errors in the various Wn and x, here denoted by Pe 
and Py: The treatment to be followed here is that given 


by Jurkevich (1976). The differential equations of the 


problem are: 


u @)\ = oe tare Ecos pat. ue (0) = 0 
see). Ou, xO) = 1 20) 
py, (t) a ee p,, (0) = 1g 
pet) 5= iecyianne te Cos aa 
m x 


where t is the phase, and gq serves to connect the stochas- 
tic £unetion wath the function being estimated. The trans= 
fer to stochastic variables 1s very simple: stochastic 


Variables will be denoted by a variable with a hat, 1.¢6., 
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ie If the differential equations above (Cinese ne. 3 
stochastic form) are integrated between observations with 
phases Cte and t, =a, the following recursion relations 


are obtained: 


a i a ) 
cae pee =, [sin*™b - sina] 
See je oe Tec cee ae! (4.32a) 
Pm(i+1) ~ Pm(iy * Paci) Osin’"b - sina] + 
FUE speo cee ay a) ee | 
Pxtl [Pe caet) ah 
Sh y(aissshy ya Sauedgenhy 2a cee ee 
f 4 py ae ; r (An 22) 
eit eit oer rears 
Mee ~ ey “oe [z,-x,] 


Here, r is the mean square value of the observational 
errors. “Equations (42324) "senve to "propagate" the solu 
tion from one observation to the next, while equations 
(45325), “apdate® the solution by processang the latest 
data point. Jurkevich (1981) has published a FORTRAN 
program which uses the Kalman filter algorithm to compute 
A £On t= 1273746 8 lhe program weseli as veryessimple, 


2m 
With eal l the “computations being. contained gin Jone DO-loop. 
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Tie presente atichor Nas adapted chs program for use on the 
MTS FORTRAN compiler at the University of Alherta. The only 
modifications required were the replacement of certain in- 
put and- output Lines. The program is listed ain Appendix 3 
The Kalman filter has its advantages and disadvan- 
tages. AS was mentioned earlier, the algorithm is compact 
and simple, and it has the desirable quality that the 
values of the AL. S obtained are least-squares estimates. 
The relative performance of the Kalman filter will be con- 
Sidered in a later chapter when individual eclipsing bin- 
aries are discussed. At this point, it suffices to say 
that the Kalman filter produces satisfactory values for 
the Ann 5° A major disadvantage is that the error analysis 
produces spuriously low numbers. Jurkevich (1976, 1981) 
has investigated this, but the cause of the problem 
remains a mystery. A possible solution, suggested by 
Jurkevich (1980), is that the system of differential equa- 


Eions, (4-32), which constitute a medel of an eEclipsing 


binary, are too simple for the purpose of error analysis. 


2m = 


AG. GBeremaAnalysis forthe A 


After the A, |S have been found by one of the pro- 


cedures just outlined, one can peniorm an) error analysis 


om S. 


LOugLVe the QiCertaintcies Ane thesvaluessoc, the A 
This error analysis may be done quite easily when least- 


squares Fourier analysis is used (i.e., equations (4.17) 


and (4.21)), since the least-squares solution produces the 
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uncertainties of the Fourier coefficients directly. How- 
even, Greer analysis 2S virtually impossible when numerical 
integration is used, since the error formulae for numerical 
integration rules are usually quite cumbersome to use (see, 
LOriexanple, Gerald, pe 2). 

DemiLcane(Cl9e lps. 127 )* nasstound= anwapproximace 


expression for the error in the Ayn 5 found by) numerical 


integration: 
AA = A ae 4533 
oY ame WS in 1 ; Ge Beil 
where Au is the error in 2 at 8 = T/2 (quadrature), and 84 
ts Vthe phase angle of first (or dast) comtact.. In most 


cases, one is interested in only an approximate value for 


AALS so equation (4.23) is sufficient for most purposes. 


To determine dA, when least-squares Fourier analysis is 
used is quite straightforward, one simply applies the rules 
Of —Errom DEOpagation to, equations, (4.4.20) 2) This, results ain 


the following expressions: 


a 
AA, = 5 Aa, ae Aa, + Aa. + 
and 
n (m) 
nN bBo ae (4.34) 
n=0 
where 
eee Ses sie Oates | 
Al) Be 2m)c. n a 1 
5 am 520 Tim+j4+a)P(m-G+5) (2541) 0,*-[nr]* 


It should be noted that the unsimplified general expression 


LOL Aon in terms of the Fourier coefficients has been used 
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in equation (4.33). If one were to use equation (4.20) 
directly, a typical result would be: 


sin“, cos“6, 
Aa + 


AA, = au Aa 
2 yg lt 2 2 2vtl1 


1-[vn/e,]° ie [ (2v+1) 1/28, 1° 


(4.35) 
with AA, s AAG being foundyin a similar fashion. «Clearly, 
only the errors in the Fourier coefficients are considered 
in equations (4.34) and (4.35). 

Considering these results, there is a definite 
advantage to using least-squares Fourier analysis to deter- 


Mine the A rather than numerical integration. The problem 


2m 
now is to relate the moments of the light curve (A,_) EC 
the geometric elements rir tos Deli peice for tne 


2 i Qua 


cases of both spherical and non-spherical stars. 


4.7 Computing the Elements 
The problem now is to relate the moments of the 


Light curve, “che Ayn (m=0,1,2,3), to the elements rp 


ay Uf 
<a 


Ss 
2m 


L,, and i. To explore the relationship between the A 


and the elements rit angen, the relatively simple 


ae 
Case On a total eclipse Of “a unsiormly bright stan will be 
considered in detail. This will serve to introduce the 
more complicated cases in which limb darkening and distor- 
tion are present. in all cases, the method of solution 
will be outlined. 


To relate the moments of the light curve to the 


geometric elements rit, ,L,, and 1 in the case of the 
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tonal Feclipse: Of asunttormiy bright. star, tolwowing the 
treatment given by Kopal (1979, pp. 148-153), one starts 


with the equation defining the Aan funcelon: 


8 
t 2m 

Aon = | Cl= Ud (Sire 0.) (42.36) 
0 


where, Once again, 95 is the phase angle of fourth contact. 


1 


In the case being considered here, one may write 1-2 in 


terms of Li and a(k,p), which were defined in the chapter 


dealing with the Russell model. The equation is: 


1-2(8@) = a(k,p) L Ce) 


7 = 


2m 


The geometric relation must also be used to replace sin” 6: 


5? = Pine he 4 + ane i " G4 3383) 


Lf 5 = cos i, (6, ato = 0) then solving ‘for sin“6 gives: 
rae = 
m —m 
eCity, 2 2 Sige 
Sine. oe mo ox) (a on) 


Differentiating with respect to a gives 


m-1 —m 
Bice) SWS Se aie Sy eS 
eG 2m . 
Therefore, Aan becomes (with (eo) = OSC ut) 
$2 ; 
m= 
eet Se ak ee ene ce (4.39) 
2m ai LD O 
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TiiLseecavwatl longus sour Starting. polnG, Lor hie snot of a 
Lotalveciipse Of a uniformly bright star. Thee limits 

55 and 65 earG the values on Oat £1rst Contact (correéspond— 
eae ise: 8.) andyat second ‘contace (j—16)). “PAtloiecstucontact, 
BiG twOestars will fappear eto ibe touching, sso 6, =r, +%,- 
Me second contact, Che larger star (with radius ro) We 
have just "covered" the smaller one (radius ri), SO ie 
met. Thererore, nT Soe oe Bor a total eclipse.) A 


May De wsplLe wpein the rollowing wey: 


2m 


63 2 
) 
om 2 ,2,m-1 .9 fl 9 gmt 2 
ra =e CSc 1 (0 = 078) ads” + (6° -67) a dé oe 
2m i 2 e) 5 @) 
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(4.40) 


The first, integral is the contribution from the total phase 
of the eclipse, in which the smaller star is completely 


covered up, with the second integral describing the contri- 


butLons £rom the partial eclipse phases. The first integral 


is quite easy if one remembers that a=1 during the total 


phase: 
® 
6 : As 
ig a m 
aces Cet weet emer Geen 
Zz O 2 
§ § 
O © 
Sete 
=m (85 om : 


Matters are not so simple during the partial phases. Doing 


the second integral by parts gives: 
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if an @xpression for 00/d6 can be found, the integration 
will be complete. For a uniformly bright disk, @ has the 


£OolLowang. form: 


2 ee ale Pang geilebe 
Us Cheer [¢, 5 sin 2¢,]+r5[¢, 5 Sin 26,] (4.42) 
Pt rt x3 Stans ees 
cos 4 = 28r, and cos Te (4.43) 


Using these equations, the following expression for 3a/364 


results: 
Ne? 2 2 2 
Y = = 
ee} esr if ay, pa ae oF: 
06 Tr, ey ‘ bane 


This equation may be written in a simpler form if one intro- 


duces a new variable 9 defined by 


Dane: 8 2 
6" = ri 2x, r,cos > + roo 
so that 30/936 becomes 
da _ _ 2 (£2) sing (4.44) 
36 T “Cy 6 ; , 
Assembling all of these results allows Am to be written in 
the form 
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result for 00/36, one can generate equations for the A._'s. 


2m 
Formmm—=.l72,0, one has (with L, =1-) 
Aa e aC | 
2 L3 
uA 2 2 
A, = Lj (c, + C5) » (4. 46a) 
= 3) 2 D) 
Ag = L, (C3 + 3050. + CiC5) 1 J 
where 
ae oan ) 
FT rise a 
Ge = ar x Bee 4 ‘ (4-54 65) 
2 Le 
Sal Ze 25 
C. = r,csc cE (eleheiis ka | 
: ; 2 2 Aces 
Solving this system of equations for rps and sin i gives 
the following result: 
B) 
S Cc 
ee (oo, and sin 4 =——_—+ (4.47) 
’ = = 
(1 C,)C, +C, Gir C3)C, +C, 


Thus, we have expressions for the geometric elements of an 
eclipsing binary in closed; form... «However, the conditions 
for suchma solution are that the disks of the stars appear 
Untiormiy sbuuaht, (net wloe realistic), and that the seclapse 
be total. The results just obtained can be made to apply 
£0 an annular eclipse if ©. and ©, are interchanged, and if 


i Z 


2 
Li is replaced by (r,/r,) Lj: The solutron for “a partial 
eclipse is an iterative one, and will not be discussed here. 
THe orocedure to be wsed 1s cuLlined by Kopal (197/97) pc .155)- 


The method of deriving the various results in this section 
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is indicative of the method used by Kopal in the more 
general Cases or aubucrary (limo. darkening and@non-spnerical 


SyestieS) ¢ 


4.8 Total and Annular Eclipses of Limb-Darkened Stars 
In MOSe Cases, the stars cOonsticlleing ai eclipsing 
binary system have a non-zero limb darkening. To modify 
Ehevrnesultes Obtained in the unitormivy braght case, to the 
limb darkened case amounts to little more than redefining 


the constants Cis Coy and, GC The method by which this 


3° 
EranslLObmation is achieved, and the resulting method of 
solution for the elements will be outlined. 


To pass to the case of non-zero limb darkening, a 


generalized form for the A. "s is required. Without going 


2m 


Into the detalls Of Ene derivation tne EOrm or the 2... Ss 


2m 
adegived by Kopal =(1979>)eq. 160) as 
Oe escon f 
: 2 ee Tico 
Som eee eT Ae DY Pvt)? (v+m+l) 
ne 2 
VaR ing 2 lr (v+tn+1L) 
x Le pe ita Gag 41,8) | 
Gala oeen errop ace) (4.48) 
n O 


where I(x) is the gamma-function, ce is the Jacobi poly- 
nomial of order n (defined as a special case of the hyper- 


geometric series oF (arb,crx), see Mathews and Walker, 
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Ogee a eC the constant defined in equation (4.29), 
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partial) follow from the equation above. 
The transtormation to an arbitrary limb darkening 
Cat De (Carried out with the aid of equations (4.29) and the 


EOL owing equations, (Kopal, 19797 pg. Lo4) 


= ee cf) 
o x Se be a v 
2S Ce) 
dq 36, os Saas (4.49) 
Gye, = C405 ue is 
where Ci, Cos C, ace Che iconstants Untroduced: an, therdis= 


GusSsion, on) the tmitormiy bright case.» In-most icasesvor 
practical interest, only a linear limb darkening coefficient 


is available. In this event, the above equations reduce to 


= ‘| 
es Gs 
15-7u 
2 2 i 
a a (4.50) 
2 2 SSS aa) | 
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These equations provide a means for introducing Limb) dar— 
kening into the method of analysis. Let us now consider 
thepanalycis Ob soclipsing Difiazy sStarsswhth)arbitcary lamb 
darkening. The Stars constitutingwthe binary (system will 
still be regarded as being spherical. 

Once, again, £ne Simplest fcase as sehat Of the coral 


eclipse. The equations required are similar to those used 
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in thevunicormiy bright case, (Kopai, 1979), og. 164 )«: 


Ay = Ly = ]-) \ 
BS) ome 
PN gi Gees arabe ory 
4 pe 2 
7 ae We So 
Ag = by (C3 +3650, + 6,85). | 


The solution proceeds in exactly the same way as in the 
uniformly bright case, once the limb darkening has been 
introduced using the equations presented earlier. However, 


the problem is more complicated in the case of annular and 


partial eclipses. The method to be presented here comes 
from a later work by Kopal (1982a). The equations relating 
Ene moments Of the Vight curve, EeEhe Aon 8! to the elements 
Pirtos and a are (for both annular and partial eclipses) : 
a | 
zo rhe 
An = i [c + (l-a Wee ear B sere 8 
eS O ee” z 
r (4552) 
ae =2, =2 oo ae 
A= L,[c, + ie (1-a.) cot i] +18, eG th 
= 3 =2- - -2 6k: Gn 
A. =1L,[C,+3c5¢,+¢,¢,+(1 a) cot | L,B,ecot i 
where Ch is the value of a at mid-eclipse, and the Bom are 


given by: 
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where me is the Jacobi polynomial and 


oye 


P- 3 
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KOpal pOintseouUte (1932, pgs 136) "that the termes containing 


the Bom are so numerically small that they may be ignored 


altogether or treated only as small perturbations. As an 


example in the same paper, Kopal computed B Bar and B 


pay 6 
for the light curve of Algol. The largest term was Bacon 
having a magnitude of approximately 5 x Tae B Cot gd and 
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B aout are smaller than B cot i by Lactors wot 100 <and) 10000 
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respectively. The method of solution to be used is depen- 
dent upon the value of aoe In the case of an annular 


eclipse, a =k*y, where k=4r,/r and where 


Dak 
gee : es a. (annular So Pee 
' a (total eclipse) 
Qo (e) 
serves to define Y. To solve for the elements in the case 


of an annular eclipse, one proceeds in the following way. 
We begin by finding ke from 


ee i 
D b 
Po Saas (4.55) 
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where ae te at maximum total eclipse and ae 2 at maximum 


tCeansit eclipse. A good) startang value tom Yas 1 Ge may 


then ber found with” the. a1q of 


a = a a 2 
A, 2p hy Bio, = Lk Vee (APx5 6) 
2 2 ce 
One then solves for x=41,csc 1 and, yY = cot: 1) from 
2 om Z ) 
_ ASAG 2A; eto ey De (te ae) ory 
= 2 2 D q 2 D 
A.A, ~ A> AL (1+3£,k°+£,k°)-3A LL, (1+£,k") +21) J 
(457) 
ell, os JN 
= di 2 
and y = A 
fe) 
iInethe equation. for x7 f., and fy are 
eee aly eRe ey 
Ean Seer ee i es 3 any) Seong as: 


These equations, along with the fact that r = kr, , provide 
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a preliminary set of elements. The value of k is recalcu- 


lated, hence a new value of Y, and the procedure is 


repeated. Only two or three iterations are required to 
provide a final set of "good" elements. The Bo Ss may be 
included if one so desires. However, it should be remem- 


bered that uncertainties in the elements and the observa- 


€1Ons May not justify the use of the Bom So Fora variety 


Of reasons, a solution tor the elements could be poorly 


determined, so the use of the B.'s would only compound the 


2m 
problem. 
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4.9 Partial Eclipses of Limb-Darkened Stars 

Ther solution for a pareial eolipsesas-cimilary to 
thats fer an annular-eclipse. In the case of a partial 
eclipse, One usually starts with k and Y both equal to 1. 


The value of a May ber found) from erther or tne following 


equations: 
= 
a6 elma ae te 5 2 
eo ME 
(42 59) 
shy Seite eee tee See 
fe) b ‘ 
Since eit Ole Ka Wu ds, =a'=a'. The value of L,. is found 
e) e) re) if 
from A, = 1 o>: One now evaluates x and y as in the annular 
Case. fo, Eind, a new value of k, one can solve for k £rom 
Rye Sia bee ay ae eo ee SE (4.60) 
4 Lis 2 fe) ' 


and then repeat the whole procedure until k no longer 
changes significantly from one iteration to the next. 

Hence we have a straightforward iterative procedure 
for computing the elements of an eclipsing binary with limb 
darkening». tor botheannularcand partaal weclipses st should 
be, noted here, that Eherdefianition- of i the tatico, of the 
radii, depends on the type of eclipse. For a transit 
ec pse, Ki-wl,. whiten fOr an woccuLballony as Je Th Ss. ainicos 
from the fact that k is defined as the radius of the 
eclipsed) star divaded by. the smadius of the eclapsingestar. 
Such, a, derination, 12s quite usefull in distinguishing eclipse 


byDeSsmi ms ties Case CL aupaltially—ecla posing systema. Lites 
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computer programs for determining the elements of an eclip- 
Sing binary for each of the three cases discussed above may 


be found in Appendix 3. 


a-e0™ Encorporatcing tne Ertects Of NOn—Sphericity 

ie che =stars Constrtucing an. eclipsing binary are 
non-spherical, one may still use the method of analysis 
just discussed, but the Aom § will have to be found with 
eherald CL equations (4721)—(4226)—or equations (4,27/)9=(431), 
i.e., the effects due to non-sphericity, reflection (and 
possibly mass transfer) will have to be filtered out before 
the solution for the elements can proceed. Other than this 
extra step in the computation, there are no added difficul- 
ties in solving for the elements of a non-spherical eclips- 


ing binary. However, one should keep in mind that the radii 


SO Obtained are in fact mean rcadiri: 


fall” Conclusizens 

This concludes the presentation of Kopal's frequency 
domain technique. The method obviously has several advan- 
tages, not the least of which is the ease with which it may 
be automated. The entire solution can be done with a 
digital computer, eliminating the need for cumbersome tables 
which are required in Russell's method. Kopal's method is 
alsomquilter compact im ts formularion, making Tia easily 
eomprehensible. Lt also has the advantage of being a 


moderne method of analysis, Since lt incorporates all current 
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knowledge of stellar structure. Moreover, all relevant 
astrophysical information can be brought to bear in, those 
cases in which one is confronted with an eclipsing binary 
having certain peculiarities. However, no method of light 
curve analysis is without its drawbacks, and Kopal's method 
is no exception. The most apparent of these is the need to 
assume a value for the limb darkening (uy) before proceed- 
ing with the analysis. One may infer a value for the limb 
darkening by considering the spectral type of each star. 
Unfortunately, the Jambedarkening inevitably turns out to 
be the most poorly determined element, no matter what 
method of analysis is used. A second drawback is the pre- 


sence of the B._'s, whose purpose, it seems, is to make 


2m 
LEee ditihicule tor the astronomer, Tndeed, one should only 
consider using the Bo 8 1f the ight curve as ob very argh 
quality (small error 22), Or 12 the solution fon the ele— 
ments is very well determined. In any other situation, the 
Bon would not be worth using. Putting these criticisms 


aside, Kopal's frequency domain technique is indeed a step 


forward, 11 Light curve analysis. 
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CHAPTER 3 


WOOD'S MODEL AND THE WINK PROGRAM 


Da De roduc eon 


tne last method of light curve analysis to be ecen= 
Sidered 15 an Gxample Of a “Synimesits Method.  Berore 
exploring Wood's model and the WINK program in detail, let 
us consider some of the general features of all synthesis 
MeLHoOdS. FUL Very simply, a SYMEnesis method is a pattern 
becOgrit2oOn Algeratinm, 22: keeps On, Constructing Light 
curves and systematically varying the model parameters 
UME Bt echieyves a match with the Gbserved Laght curve. 
The algorithm synthesizes the light curve using a model of 
the eclipsing binary, hence the name. Computer programs 
sing the Light curve Synthesis approach are wsually large 
(eyprcally 2000 lines of code) and GUuLte complicated. Most 
O£ the currently available programs are the end products of 
several years of development. 

From a more Eechnical pormt Of view, Synthesis: pro- 
Grams Wellazge either a triaxial elilaosoid cr Rocke suxrtace 
Geomecry fOr Ehe Component stars. Wood °s model wees tie 
former. A typical synthesis program proceeds in two general 
Sreos) Ene caleulation of Ehe Lumimosity 2 at all observed 
Phases 6, and a parameter adjustment step to bring the 
model into better agreement with the observed Light curve. 
Ths step 1s Ustally achieved wht the atid ot a Girtterential 


Gorrector operating im roughly the Same way as in Russell's 
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method. The stellar models are usually complex, incor- 
porating the properties of the atmospheres of both stars to 
determine parameters such as surface temperature, surface 
Gravity, and limb darkening. Leis also possible to accouns 
for extended atmospheres, which are found in Wolf-Rayet 
binaries and systems that have red giant components (a 

well known example of the latter is VV Cephei). To handle 
particularly close eclipsing systems, a typical synthesis 
program would incorporate a subroutine to handle the reflec- 
tion effect. There is no standard procedure for handling 
reflection, and the approach used in each computer program 
is different. When dealing with the Pe gieecion effect, one 
has to make a tradeoff between computing time and the 
accuracy of the reflection model. Reflection effect cal- 
culations can Significantly increase the running time of a 
synthesis program. Clearly, synthesis methods make use of 
avcomputer ss number=crunching? ability... This allows, one to 
use the complex stellar models mentioned earlier. The first 


computer programs using the synthesis approach appeared in 
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the early 1970's. Every researcher had a different approach. 


The earliest synthesis programs were those devised by Wilson 
and Devinney (1971), Hill and Hutchings (1970)" and Gucy 
(1968). In each case, the approach taken was approximately 
the same, namely to use a Roche geometry, and to compute 
the surface intensity I(u) at several thousand points on the 
surface of each star. These programs were closely followed 


Byruhioseson Wood (19727) Rucinski (1973), Mochnacks sand 
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Doughty (i972), Nelson and Davis (1972), and Berunier (1975): 
Mochnacki and Doughty were able to simplify the synthesis 
procedure’ by using cylindrical coordinates rather than the 
usual spherical bipolar coordinates. Berthier's program 
used a library of over 4000 known light curves to find ini- 
tial values for model parameters. The Nelson-Davis approach 
uses an unorthodox integration procedure for determining 
surface brightness. More recently, Hill (1979), Binnendijk 
(1977), and Budding and Najim (1980) have proposed synthesis 
methods, with various improvements such as faster integra- 
tion procedures, and an increased use of analytical formulae 
for computing certain quantities. Binnendijk's program even 
allows the user to compute a theoretical radial velocity 
curve for use with spectroscopic data. A paper by Hutchings 
(1971) describes the computational aspects of a synthesis 
method in some detail. Hutchings also discusses the compu- 
tation of line absorption and emission profiles for eclips- 
ing binaries, as well as the treatment of extended or 
expanding stellar atmospheres. The particular model that 
Will be discussed here is Wood"s model,in part because it 1s 
one of the most widely used programs. The WINK and WINK8 


computer programs, both devised by Wood, will be discussed. 


5.2 Model Parameters 
We begin the survey of Wood's model and the WINK 
program by describing the general properties of the model. 


Tie treatment Given nere L.Ollows thac~n Woods (197 171072). 
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Tie node l weakes sini account “Otarional and, eldaal  discor— 
tion, limb darkening, gravity brightening, and reflection. 
The three major simplifications employed are that the stars 
are etriaxial ellipsoids, with the reflection effiece being 
approximated, rather than being computed by more rigorous 
(and consequently more time-consuming) means, and that the 
stars rotate in their orbital plane with a period equal to 
the orbital period. The first of these assumptions provides 
for fast Computation, and is valid except for very close 
systems since the Roche surface defines the extent of the 
outer atmosphere and not the photosphere. Wood claims that 
the approximation used to simulate the reflection effect 
(to be discussed later) is "very good". This is to be 
interpreted in the sense that the approximation agrees well 
with more rigorous computations (i.e., the work of Napier 
(1968)).. “From an astrophysical point of view, the rerlec— 
tion approximation is a grey atmosphere approximation. The 
assumption regarding rotation can be relaxed somewhat, since 
orbital skew and polar tilt are allowed. Orbital skew 
implies that the stars do not face one another along their 
Mayor axes, but that the major axis = of one Of the stars 
either lags or leads periastron passage by an amount o. 
Polar tilt allows for Stars whose rotation axis 1s not 
parallel to, the orbital pole (by an amount 1) Inveach 
case, the deviations must be constant. 

The model is defined by three sets of parameters, 


namely the orbital, geometric, and photometric parameters. 
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The orbital parameters serve to specify the nature of the 
relative, OLbLe Gf Star Boabout Star A. Star A is eclipsed 
ducing the primary eclipse (i.e., it is the primary compo— 
nent). More specifically, the orbital parameters are the 
period P, time of conjunction ER (usually zero), the semi- 
MajOpeaxts OL the sorbic Ry (usually one), eccentricity e, 
fongutude Of periastron ©, and inclination i (see fig. 14)": 


The definition of the inclination remains unchanged from 


previous chapters. The geometric parameters are just the 
stellar semiaxes a, b, and c. Certain relationships that 
exist between the semiaxes will be explored later. The 


photometric parameters specify the apparent intensity dis- 
tribution on each star. These parameters are the surface 


intensity I, or the value of I, ule’ 
I= I, (1 = a + wi cos) 7) 


at time opto ees (time of quadrature), the limb darken- 
ing CcOecbiicuent UW, Gravicey brightening cocriicient vy “and 
reflection coefficient (or albedo) w. 

As mentioned earlier, there are certain relation- 


ships between the stellar semiaxes. Wood replaces the semi- 


axes a, b, and c by a new set of quantities 


an = aR, 
an = k aR ? Peers, ke = an/an 
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FIGURE 14, ORBITAL PARAMETERS IN WOOD'S MODEL, 
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The parameter ne “US SO deraned that, Le aay then the pri- 
Many eclipse 1s) an occultation, and a transit, at ee 1. The 
Stars aye assumed to be triaxial ellipsoids. Chandrasekhar 
(1933) showed that the stellar semiaxes can be written in 


terms of the polytropic index n, the mass ratio q (mass of 


Stam e/mass Of Star A) and a quantity \ = a/Ror where ay 
is the unperturbed radius" of the star, i.6., the radius 
Olea sphere having a volume <qual to that of the stars If 


each of the semiaxes is expanded in a series in v, n, and 
q, one obtains the triaxial ellipsoid geometry (retaining 
CeErms up cOVorder 3°in Vv). 2p 1S cenvenient to antroduce a 


quantity ue which is defined as 


One can now express either the semiaxes or the quantities 
iWeroduced in equation (Seek), in terms tof Aca! Ka q and nn 
As in the models discussed in previous chapters, the 


intensity at any point on the apparent disk of either star 


may be written as 
I = T(t = Ue COSeajin, (CSB) 


where I, is the intensity at the “sub-earth” point, or the 
POLnt at whieh the line Of ssight ais pespendieularn co: tne 


stellar surface (y= 0). in the original WINK program 
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(version published by Wood in 1972), Wood uses a blackbody 
approximation for Io: Improvements on this will be dis- 
eussedslatera The Gravity bGrghtening v. si defined, in she 
following way (analogous to the definition of limb darken- 
ing) 


Dao eee es GMa eal (5.3) 


where r is the local radius (at any point), and r is the 


radius at the sub-earth point at time T I was defined 


0° 
earlier as a photometric parameter. 
The quantity which one must compare with the obser- 


vational data is the system luminosity at a given time t. 


To compute this, one must evaluate 


ity [faj+r) a- Bh a2 bh elo Wel 4, (eyrare) 


the integration being taken over the apparent ellipse of 
either star (shape of star projected onto the plane of the 


sky). The system luminosity at any time t is 
L Ce Me Re ibs Wey) = Soi (ie), (55.5) 


is the light Loss during eclipse. Ly and L, are 


= Se =1. The WINK = 
Q! oe La L, tk e DEO 


Gram allows for the presence of a third star (the =third 


where L 
ecl 


SOxsderinea that at t= Er 


ligne) in the Lollowing sense: 


F = 1 
L, + L, La 


Hence, the calculation of the light cumve amounts to the 


computation of Leo ft): In general, three integrations are 
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required, one over each star, and one over the overlapping 
area gr any). «One must, know the outline. of each stam and 
the overlapped area,..as well as the antensity at any point 
Om wene. 6itars.. Also, reflection must. be: taken. into account. 
The contribution to I from the reflection effect is the 
quantity i> that appears in equaizr1on (5.4). Wood computes 
ny ay determining. Ehe, local incident jatensity &” at a, maiven 
point and then reflecting a fraction w uniformly over the 
outgoing hemisphere: 

ic, 


rt = _ wi mre (5.6) 
2h wel = 5) 


i mist Wow be added “CoO qa: The method used to compute Le 
wil beldiscussed Gn detail later; but the approximation 

mwsed 25 valid for an and an less then 0.5 (fhe usual range 
Ge interest). Weed also eailows for the possibility Of. an 


extended atmosphere around either star, but this aspect of 


the model will not be discussed here. 


3 Computation of the Blements 


Me solution for the elements rs cone by least— 


Squares differential corrections. The hight curve iIntensuty 


1 2S assumed £o be of the form T=I(t, Xr -227 XO), where 
Xs are the elements (unknowns). We may then write 
a1 a1 
= = —— AX +..4 ——— 
of ee cade aX ak aX AX 
a N 
N 
aL 
ent Ar Se 
e) ok Ps f ( ) 
i=l A 
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where AX. iS the differential correction to the estimate 
X.- Fach ODServabion provides an equation of condition 
having Ene £Orm Of equatven (5.7)... The resulting System of 
equatrtons 15 solved by the least-squares method. This may 
be turned 1ntO ah Puerative process whach, one hopes; will 
converge =O a solution. Convergence 1S net easy to attain 
Since there are observational errors and @gnterrelationshios 
beeveen Lhe parameters. It is @ieo GLiticult to know whieh 
pPawameters should be varied, and a Lot of guess work as 
involved. In géneral, more than one computer run is re- 
quired. The partial derivatives aL/oXx, are computed numer- 
ically, and there exist certain short-cuts and simplifica- 


EROS 2n Lely Computation. 


Det Detalles of the WINK Program 


We now consider some specific aspects of Wood's 
model and the WINK program. Of greatest interest are the 
Merhoa Of Aanteqration used, the way in which eclipses are 
detected and thein limits* found, - the wet lection approxima-— 


tion, and the use of model stellar atmospheres. 
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ity Order to discuss the integration precedure used 
in the WINK program, we must first define the coordinate 
systems in which the dmbtégration Ws Carried oul. "The fun- 
damental coordinate system is a rectangular one centered on 


Star A, The s-axis points along the line of signv wand the 
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V2 plane 2s in the plane of the sky. Hach star has 1ts own 
coordinate system, denoted by We"py',z'). The Z"=axis is 
coincident with the star's principle axis. Since there is 
octant symmetry, the sense of the coordinates does not 
matter (except for reflection effect calculations). The 
apparent ellipsoids each have a yz coordinate system lying 
along the major axes of the apparent ellipsoids. The yz 
system is rotated with respect to the yz plane. 

The various integration procedures used in WINK are 
handled by the subroutines TOTINT (total eclipse), ANNECL 
(annular eclipse), ECLINT (partial eclipse) and ATMECL 
(atmospheric eclipse; replaces TOTINT). Subroutine TOTINT 
alse computes: they totals lighteoutpuu frompeach Stansac any 
time. In all integration computations, Gaussian quadrature 
is used, with three grid sizes available: 4x 4 (coarse), 
Geo i(normaleor- detault), ands 12x12. (hugh) precision) ..),00 
integrate over an ellipse with semiaxes a and b, one uses 
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» W. we (ax, ; X, bv 1 = a) 

where Bee) is the intensity along the line of sight at 
pointe \y72)) (equation (5.2) )=, W. and ie are the Gaussian 
weights, and Xo = are the Gaussian ordinates (loaded by 
Subroutine GRID). Also, m 2s the number of quadrature 
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a we = a 
ee ae ey ‘ ee a 27, 
1 a 
ea ops Care 7 De eae ce: 
Yu and Y, ane the) Vy-limits oreanteqration (Y., - Ys found 
by subroutine LIMITY. Zr and Za. are the z-limits as found 


by subroutine LIMITZ, and are functions of ae 

In all cases, the integration is performed in the yz coor- 
dinates, except for ECLINT, which uses a coordinate rotation. 
A coordinate translation is required in ANNECL because the 
integration is performed over the area of one star with the 
intensity poings of the other star. Acvsimilar situation 
exists when star B is eclipsed, since the integration limits 
are in the coordinate system of star A. The intensity 
Laly,2)) 1S computed bystene function BRIGHTS As andicatced, 


p 
BRIGHT uses the yz coordinates. Figure 15 shows the various 


ineeo vation .rids. 


5.0 HClipse Detectionvand Limit Handing 

hes thesduty of ‘the subroutines, SCREEN, LIMITY, and 
LiMEUZ | tO.search for and find the lamats of any eclipses 
thatweccur during computation of the light curve Subroutine 
SCREEN searches for eclipses by comparing the centre-to- 


centre separation $6 with the sum of the radius vectors. 
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Obvicusiy, “an eclipse occurs if 6 < rj t+r,- Subroutines 
LIMITY and LIMITZ determine the y and z limits of the over- 
lapping area during an eclipse and the nature of the eclipse 
(Lies eOcalpannilar, Or® parneial) In general, two 
ellipses can intersect at.as many as four points.  LIMITY 
(and LIMITZ) uses a bisection search (a line y=constant) 

to find these intersections. For an eclipse to occur, each 
ellipse must, have two. rocts (four roots in total). The way 
in which the roots interleave determines whether an eclipse 
SCCUES, and 2 it is partial. The roots are denoted by Z 


Al 


and Zn5 EOr estar 7A, and #by 251 and Z 25 fOr Star Ba yan 


eclipse does not occur if Zn and Zn are both greater than, 


or both less than 251 and 239° Figure 16 shows the various 
cases that can occur if an eclipse is suspected. The bi- 
section search (Wood calls this a "Scan wire") proceeds 
from right to left to determine where the interleaving of 
the roots changes. When such a change is found, the search 
is reversed and a smaller search interval is used. If the 
procedure is repeated, the location of the intersection may 
be determined quite accurately. In the event of a shallow 
eclipse, a smaller search interval is used, and if such an 
eclipse occurs On’ Enel negative y-axis,)the negative Aims ¢ 
PSsetound Lircst, Or in Other words, the’ search goes Crom 
leper co rights. The bisection algoricim tsedemenseraved 2m 
Ptogmpenl /saeA) Salil at i DCOceOuLecmissUseds LOnUeLeLM nem rie. 2— 


(iis ronmehereoclLipsce,., and it shoul dsbe noted tiat thes 2— 


limites are usually the two. "inners 1oecs (see, figure 16) 
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2 ROOTS - NO ECLIPSE UPROOTS = NO ECLIPSE: 


4 ROOTS - TOTAL OR 4 ROOTS - PARTIAL ECLIPSE, 

ANNULAR ECLIPSE STAR LIMBS AS INTEGRATION 
MIS: 

4 ROOTS - PARTIAL ECLIPSE, 4 ROOTS - PARTIAL ECLIPSE, 

INTERSECTIONS AS INTEGRATION GENERAL CASE OF FOUR 

Mss INTERSECTIONS. 


FIGURE 16, IDENTIFICATION OF ECLIPSE TYPES (ADAPTED FROM 
WOODKLO/ 2yaPGy 20). 
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5./ The Reflection Effect 


the reflection effect 1s nhanddeqd- bythe subroutine 
KREPD.= The general procedure for determining the reflected 
intensity I* was discussed earlier. Let us now examine the 
computation of L*, the incident intensity. L” is assumed to 
bea function of the ‘source star’s apparent ‘angular size 
(AL, as "seen £rom a point om the other star), the inten-— 
sity Te at the end of its a-axis (1.e., the sub-stellar 
point), Lts* Limb darkening a ic and its apparent zenith 
distance (\"). 9 "The approximation used was derived from 
more exact results obtained by numerical integration, since 
such an integration can increase the computing time required 


by an order of magnitude. “We may now write L* as 
d £ 
L = TLA, (ul) g(cos i ae 


The limb darkening dependence, which appears through fF(ul), 
is well approximated by 


aS 
Se pnb ae 


Since 


(l=u+u cos yy) siny dy = 1 -— 5 
0 


In other words, one is considering the contribution from 
the visible hemisphere of the source star. Wood has found 


that g(cos x") can be approximated (linearly) as 


g(cos A') = -0.065354+ag, + (2.044 + a,g,) cos Ns a © 


where 
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g, = 0.224935 - 0.761696 a_ + 3.81425 at 


Go = 0.170831 - 1.231707 a. t 9.955083 at ' 


anc AG tsowthiewlarger of 0i-or 


Onc Gt a (-0.82442 ot 1.4343la_) + 


SF (cos A') [-1.22172 +a, (-0. 43316 + 4.9378 a_)] 


In the above formulae, a. and a are the apparent sizes of 
the source star and reflecting star respectively. The 
equations and procedure for determining \' may be found 


in ythetpublication by Wood (1972, Appendix 2) . 


5.8 Model Stellar Atmospheres 

The original WINK program, as published by Wood in 
1972, uses Planck's blackbody law to simulate the atmos- 
phere of each star in an eclipsing binary system. In 
general, stars are not blackbody radiators, but over a 
narrow wavelength region, the approximation is valid. To 
allow U, B, and V light curves to be analyzed, one must 
have a better model. of a star*s, atmosphere. To this end, 
fluxes computed from LTE (local thermodynamic equilibrium) 
model atmospheres were incorporated in succeeding versions 
On the WINK program. The current version, WINKS; uses the 
LTE model atmospheres computed by Kurucz (1979). By speci- 
fying the wavelength of observation (in Angstroms), the 
etfective: temperature of each star Min degrees, Kelvin). and 


the Logarithm-of the stellar surface gravity (g), WINKS 
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Can interpolate in a table of fluxes to produce the 
required flux (and intensity) for the given wavelength, 
eflective Lemperature, ands logec. ssThis procedure is 
handled by the function ATMOD. The atmosphere tables 
Govetva tange Of Sumtacescrava cies (log sg, = 62.45 seo 64.5 )u, 
wavelengths (3300 A to 8000 A), and two temperature ranges 
(4000 Ke to, 9500. Kvand 11000.k.to,40000 kK). In the compu 
tation of the atmospheres, Kurucz used more realistic line 
opacities to simulate a wide variety of stellar spectral 
types. The atmospheres were also computed at 342 different 
wavelengths. Clearly; an approach) suchas thas has a 


definite advantage over the blackbody approximation. 


Da) BCOnCLUS TON 

Thas beings us-to the end sof the discussionsof Wood's 
model and the WINK program. Evidently, synthesis programs 
such as WINK and WINK8 use a much more realistic model of 
eclipsing binary systems, constituting a definite improve- 
ment over the Russell model. The stars are no longer 
restricted to a spherical geometry, so that one approach 
can be applied to a wide variety of systems. However, in 
certain cases, such as the W Ursae Majoris stars, one must 
EeSOreE to the Roche geometry. In contrast, Kopal”s Ere- 
quency domain method requires two different approaches, one 
for celacively “wide binary systems, the other (tor close 
binary systems. Another advantage of a synthesis program 


is that complex or subtle effects can be handled with 
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relative ease (with the differential corrector), once the 
more fundamental parameters (such as radii and inclination) 
are known. Furthermore, the reflection effect is handled 
ina Satisfactory Manner, not only in a physical sense, but 
also ins terms of, computing, time. \Finally, realistic «stellar 
atmospheres are used, allowing one to analyze the light 
curve at different wavelengths (i.e., BV or UBV). Model 
atmospheres such as these also allow a convective atmosphere 
to be simulated. Other atmospheric peculiarities such as 
starspots or extended atmospheres may also be included. The 
only drawback worth mentioning is that convergence to a 
solution is not guaranteed, and one must beware of non- 
globals minimal (4..e., ee) may have several minima, only 
one of which is the true solution (the deepest one)). How- 
ever, good initial estimates of the elements and good obser- 
vational data can alleviate most ambiguities. Clearly, 
models of eclipsing binary systems such as Wood's model 
provide one with a more detailed and physically realistic 
model, which is something of great necessity when one is 
dealing with eclipsing binary systems having photometric 
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CHAPTER 6 


SOME PRACTICAL GUIDELINES 


FLOM “Ene preceding discussion, te is clear eta 
some of the models of eclipsing binaries can only be used 
in certain Situations, and that each requires a different 
eamoune Of Initial data. “For an inttial selution, tne 
Russell-Merrill method, or the version due to Tabachnik, 
may be used. This will provide one with a set of "rough" 
elements, which may then be used as initial data for a 
synthesis program such as WINK8. Moreover, the Russell- 
Merrill method requires only one-half of an eclipse, and 
may therefore be used if the light curve is incomplete. 
The preceding remarks also apply to Kitamura's method, since 
1€ iS a version of the Russell—Merrill method. Solutions 
made uSing either Kopal's frequency-domain technique or a 
synthesis method can be regarded as final, since both 
methods model a binary system in some detail. Kopal's 
method requires only one-half of an eclipse, but if there 
is a significant out-of-eclipse effect, this portion of the 
light curve is required. Since no initial estimates of the 
elements are required, one may make a solution directly from 
the data in one step. On the other hand, a synthesis pro- 
gram requires a complete light curve (implying the need for 
over one hundred data points) and an initial set of elements. 


However, the corresponding solution is (in most cases) of 
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high quality. Therefore, Kopal's method and the synthesis 
programs should be regarded as competitors. At this point, 
it 1s worth describing the procedure to be used if one is 
analyzing a close eclipsing binary. In this case, one 
Should avoid the Russell-Merrill method, or any of its 
VaLlancsy sineeytnevunder | yang physical modelmis nour aieail 
realistic (e€.g., the lines of constant brightnéss on the 
SUBLaACe- Of an ellipsoid are; note econcentyic, eOntrary co the 
assumption made in the Russell model). Moreover, such an 
analysis would require the use of rectification, which is 
not a valid procedure. If one is attempting to analyze the 
light curve of a close eclipsing binary, either Kopal's 
method or a synthesis method should be used. To provide 
initial data for the latter, one may use a Russell-Merrill 
method since any effects due to either stellar distortion 
or to matter streams are at a minimum during eclipse. In 
fact, WINK8 provides its own set of default elements if no 
initial ones are available. In using Kopal's method, one 
must use the procedure described for the analysis of close 
eclipsing binaries “(equations = (4521) co (426), sone 4 14) 
Sic (Aer weco 14) sol) wor enapre rea) mem DUied MG, ali dC hie selec 
avialysis;, UeMis advisable to assume a constant Values fer 
Ene lamb darkening, Since = asvariacionwin’ theslim> darkening 
Canemamuc a change Ln Ene radius Or seilther Stan (see kopal: 
(E979); pge 232). These comments should serve as -asbroad 
guideline: in the analysis of any eclipsing binary ,)buc rene 
Shoulda also be Caretul fto account 20r che wndavidual pecu— 
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CHAPTER 7 


Hoo HERCULES 


i. 2 ner oclie ti on 


ihe first example of the analysis of a “light curve 
OLvan eclipsing binary will bewtthat of, HS Hereulis. » The 
Gate tO be used in this analysas come from the paper ct 
Hall end Hubbard (1971). The results obtained in this 
chapter will be compared to those of Hall and Hubbard, who 
ised the Russell —Merrili method (wreh) recri :icatwon) ito 
obtain the geometric elements. A-light curve appears in 
froure 18. 

| The most noticeable feature Of HS Hexrculis tsAthat 
the primary (deeper) eclipse is annular, and that the 
secondary is total. An annular eclipse is identified by a 
rounded minimum, unlike the sharp minimum of a partial 
eclipse, or the flat minimum of a total eclipse. Therefore, 
the larger star is also the brighter of the two. in most 
cases, the situation is reversed, with the smaller star 
beang brighter. Since the secondary (total) eclipse is 
well-defined, it may be used in determining the geometric 
elements. Another feature that can be detected after some 
inspection is the displacement of the secondary minimum 
toward the primary minimum by an amount AG = 0.02 revolu= 
tions. One would therefore suspect a non-zero eccentricity 
and significant apsidal motion. Hall and Hubbard propose a 


value e=0.033 for the eccentricity and an apsidal period 
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P,= 15-5 yr. Scarfe and Barlow (1974) were unable to con- 
firm thissapsidalsperiodam ihe values of chand P, must be 
considered doubtful for another reason. The present author 


has discovered that the heliocentric Julian dates and 

phases as given by Hall and Hubbard do not coincide with 

one another (using Hall and Hubbard's values for P and ty): 
Theresissa constant error of \-0,0009 (+1 in the last place) 
revolutions when the phases of Hall and Hubbard are compared 
with those used in the present work. The cause of this dis- 
crepancy is not known, although the possibility of not making 
the heliocentric correction to the geocentric Julian dates 
hasewmec, ‘been ruled cut. The light. curvershown in figure 
18 incorporates the present author's phases. The two stars 
constituting the HS Herculis system are of early spectral 
type, with the larger one being of type B4 and the secondary 
of type A4, as given by Hall and Hubbard. The final feature 
worthy of note is that the secondary minimum is asymmetric, 
Since the system is brighter after secondary eclipse than 
it is before the eclipse. This effect is most pronounced in 
the U Light curve, but is barely, noticeable. am the V daght 
curve. This asymmetry has been attributed to a gas stream 
near the inner Lagrangian point of the system, which is 
Situated in such a way that it is inyisible during the des= 
cending branch of the secondary eclipse, but reappears as 


the system comes out of eclipse. 


Hhater ehionies 250 ob. Sisadul | 
«fae bis > sob epulew 2 hcaddut, chide 
lerety teat ee ne PR) C000.ee iG ns ements 
herenyace Ste Headeuit brits c5syl 0 apsadg arte nei whos 
268 2ila-Io sees SAT. tee gama al wk Seay @ 
Salter con to yitlidigeug arty igpodoks «avehs Jom t 
cegas nabiet ofasaepbey fs oF rolsodwuep oie 
euded? ul duodk avuto asit,ek? ..dve Seles need fanici 
erast ows et? .Beakiid a2! ondgce deegty O09 aps s SeRyVOUL 


bextooada ylike 30: wTs mates <cleto7e4 22 od? gat 


Ytcohacase sity fac OF SGyS To. paesed eo yepzel ony athe 
avvdbet tent of . Sxteddely begin ed devi ee abA 
chaseniyes wa: attain Piphneves Bh) eens Gt oidn bo’ ydato “ae 
nies ougéioca rR oan cI «syte teatipfat ef megeya ent eond i! 
ni toneohous teem el ssetie say -seegnice ects sroted a - 


jipet 9 of-aP didverissr ¢lssmn sl, ane -ev i> attgtt 0 od 


— 
= 7 


onade 2ny 2 OF Baivdintts ceed sac yoda ain?  evaus 
2) dordw ,mebeya ed? 16 Juccg he bene ytd yseurel eda. 3 pt" 
~est aii. obixab aldtesves a: 2: oee vow @ dove Rake 7 
ve shedgeney aud wank ba eieionie ale LO donecad ® 
sedi ine To duo Semen aX 


- 


LOG 


Leo suigGhte curve Analysis 

The geometric elements, for the methods of Russell 
(Tabachnik's method), Kitamura, and Kopal are presented in 
table 7.1. The results from WINKS are presented in table 
7.2. In computing the geometric elements with the first 
three methods, the depth equation™ (equation (2216) ) was 
SO .VeCe .Olsk, thes ratio: Of the radia. thas was) done using 
the iterative procedure suggested by Jurkevich (1970, pg. 75). 


In this procedure, one solves the depth equation for oe 


1-A 
ay b 
k = Senge iy Gia) 
a 
where 
i 2 
Y(k,-1) = [3(1- x,) + 2x, ( GNGe) Ale WS > say) 
and 
di 2 ~ cl ee 
T(k) ==—[3 sin Vk = (3 —4k) (1 + 2k) VE - kk). 


Shih 


Here hy is the value of 2(8) at the moment of internal tan- 
gency during the transit eclipse and Se is the depth of the 
eccultativon eclipse. Also, x is the limb darkening of the 
larger star. Equation (7.1) is solved using a simple iter- 
ative BELOceaure.. A prhoOgrameron ther Ti—59  oregrammnabie 
Calculator, whach performs “this calculation, may beyiound 
in Appendix i. This’ procedure: preduced) anvalue ok k= 02532 


(fOr = 0 260)) whichis sliionely larcern= chan sali and 


b 
Hubbard 7s 02 55% wn all calculatirens;.-a imo derkening 


Xp = 0.6 was used. 
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Russell Model 
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Sum of squares of residuals 
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The ‘Geometric Hlements of HS Herculas 


(Tabachnik's method - total eclipse) 


corrected Gorrections (ZAPP.DC) 
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Mable 7 sl (Cont “d) 


Kopal Model (annular eclipse) 


ry (eclipsed star) Or OR 07, 006 
vr. (eclipsing star) O26 520.024 
k NEC 2 03S 
i C2 ogee oye 
x 0.60 (fixed) 
Ly (eclipsed star) 0.902 (fixed) 
L. (eclipsing Stax) 702098) (fixed) 


Moments of the light curve (program EB.FS): 
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Jd.2 WINKG. results for HS Herculis 


a) 
Il 


O26 0003 


1 
Kees Do at 000s 
ry = OLE 4 Smee O03 
L@e=e 0072 0°RS 


Magnitude at quadrature = 8.522+ 0.002 mag. 


KF 
II 


23.000 = 600°K 


ral 
Il 


LOCOCO rR Kt used) 
x = 0.6 (£ixed) 


log 3, = log, = 4.0 (fixed) 


q = mass ratio = 0.3 (fixed) 
r2m.s. CLLOor = 20.013 
V(0-€)- = 0.0227 

Number of iterations = 4 


Astrophysical Parameters: 


Vv es Hoa! K) ote. K) Polytrope 
Oreo One2ow 2O0G 225 23133-52109 ey 0) 
0.145 Oa 10000. 00 10064.91 5.0 
a=0.266 k= 0.553 k= 0.553 Jango ioe 
J,/I, = 0-25 daira = 0.22 Jn 6 = 0.035 

ellipticity ig U,=x Vv Ww a b S 


029921) =050036 (0.76 f=4.0 0s 0 2668 02640226) 
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Table j.2 feo tad) 


LUMLHOSLELEeS : 


Star apparent normalized total (47) 
A Onis 0.929 as sato 
B Or Ones 4 OS0 Zbl deo24. 5) x ioe 

ratio O07 65 020106 

Reflection effect: 

star unheated (°K) heated (°K) 
A 22885.96 DREN ERENEY 5s SRV) 
B 9845.18 13709206 


normalized 


O98 955 


0.01045 


(47) 


LO 


ale adh 


Tabachnik's method (program LINE2) produced slightly 
smaller radii and a lower inclination. The most probable 
cause: of this is that rectification was’ not used, since out-— 
of-eclipse effects (other than the gas stream mentioned 
earlier) are minimal. One would expect rectification to 
alter ‘the sizes “of the "stars and the anclination by aiemall 
amount since rectification produces an equivalent "Spherical" 
system. Kitamura's method produced initial elements very 
close to those found by Hall and Hubbard. A differential 
correction (using program ZAPP.DC) did not produce any dras- 
tic changes in the elements. Kopal's method (programs KAL, 
ANNULAR, and EB.FS) was used with data from the annular 
eclipse. The elements found by ANNULAR are close to those 
found using Tabachnik's method. A further analysis, using 
programs EB.FS and ERROR, established the approximate errors 
in the elements. Program EB.FS (see Jurkevich (1980)) fits 
a Fourier series to the eclipse data, thus determining the 
Fourier coefficients and their uncertainties. It then 
determines the moments and their errors (see table 7.1). 
Program ERROR then uses these uncertainties, and the known 
values of the elements to compute their errors (see Kopal 
(1982a), pp. 154-155). The WINK8 solution used 143 data 
points, evenly distributed throughout the light curve of 
HS Herculis. The elements obtained were identical to those 
founds by wail fand Hubbard. > lt-iaseanteresting Veco oce chat 
only four iterations were required (a maximum of six was 


allowed)> “and that Gn the last iteration; “the anclinarion 
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wasrdeleted from the list of variables to be corrected. 
WINKGE will do sthis 4f 4 ‘cormmection AX. is lower than a 
certain preset limit. One may conclude that, in the case 
of HS Herculis, the solution was very well determined. 
Reflection effect calculations revealed a significant con- 
ErabupionsLromethetpraimary W(notcer star jira. -ethe  ttem= 
Peracure vor “the secondary (cooler) star ssw signiticantty 
affected by radiation from the primary. However, the con- 
verse is not true. However, in the WINK8 run, one temper- 
ature was held constant (T, = 10000 K) and the other allowed 
to vary. This would influence the magnitude of the reflec- 
tion effect, so the result should be treated with some 
caution. Also, the temperature of the primary obtained by 


WINK8 does not agree with the result of Hall and Hubbard 


Keotwa B4atstar LT] L6000TK). 


(/@smeeCOncLlusions 


in concluding the analysis ofl His Herculis ; chere 


are a few points that should be discussed. The most notable 


feature of the various sets of elements obtained is that 
their values are similar, but not the same (in general). 
These differences arise primarily from the different models 
and computational procedures used. In the case of 
Tabachnik's method, one is fitting a line to the variables 
Since and Gueacoy so that ‘the schoice of k hasan eittect 
on the elements obtained. Another factor influencing the 


solution is the scatter of the points in the eae 5 - (Geeks) 
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plane, since a least-squares line minimizes the square 


distance (l.e., (Grerens) between the line and the data 
points. This is important because the elements depend on 
the slope and intercept of the line. Since Kitamura's 


method involves the use of the Fourier transform of 2(6), 

a smoothing operation is involved, decreasing the influence 
of scatter of the data points. Since Kitamura's tables are 
quite comprehensive, a good first approximation to the 
elements may be obtained. With reference to both 
Tabachnik's and Kitamura's methods, one must remember that 
they are "indirect" methods, since the elements are obtained 


2 


through the use of Since Gest) Paine i and PL/F 5: 


2! 
Therefore, elements differing by small amounts should be 
expected. Furthermore, it should be kept in mind that the 
Russell model, as used in the methods just discussed, is a 
first approximation only. The model is almost entirely 
geometric in character. The only physical features used 
are the limb darkening and the relative luminosities of the 
stars constituting the binary system. The essential point 
is that a simple, approximate model will give approximate 
elements. Since Kopal's method uses an entirely different 
approach to the problem, one would expect some differences 
between results obtained with it and with Russell-type 
methods. In Kopal's method, the limb darkening enters the 
problem directly, So One has amore, piysicaily ealasiure 
model (one can easily include non-linear limb darkening) . 
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Operation and a subsequent error analysis. In complete 
contrast to the preceding methods, WINK8 fits the whole 
light curve to a model, whereas the previous methods use 
only eclipse data. Moreover, WINK8 uses a much more realis- 
tic model, in which the stars are treated as ellipsoids, not 
spheres. 1t also allows for finer details such as reflec— 
tion, atmospheres (hence temperatures), both linear and 
non-linear limb darkening, and gravity brightening. There 
are also other parameters which may either be varied or held 
constant, according to the needs of the user. From the 
results obtained with WINK8, one can see that the components 
Gb Shs Herculis are spherical toma very Good approximation. 
Therefore, when one is comparing elements obtained by 
various methods, it is important to understand the assump- 
tions and procedures involved in a particular solution 
“qevhed: 

The next point to be addressed concerns the identi- 
fication of the "best" set of elements. Technically, one 
should choose that set which minimizes ee However, 
we must remember that the methods of Tabachnik, Kitamura, 
andiykKOpal ise only eclipse data, and chau winks: “ses ehe 
whioleulight curve, as mentioned carnlien....One Ccannoustehen 
compare the other three methods to WINK8, since these 
methods assume a constant value for %(8) outside eclipses 
(in “the spherical approximation). The methods of Tabachnik 
andakleamusa appeatmetOuGlVewthioubDest& Fit) tomeiemeclipse 


data (by considering Cen as while WINK8 provides the 
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npestetit' Lo the whole light curve. “One might then ask 
which method provides a set of elements closest to the set 
obtained by WINK8. From this point of view, the methods of 
Kitamura and Kopal are the closest. It should be noted that 
in general, the inclination has a greater uncertainty with 
respect cO=the ovher “elements. "This can be seen im Vvarvous 
sets of elements obtained. Therefore, one can only choose 

a "best" set of elements after an intercomparison of the 
various sets of elements, and after considering the value 

of J (0-c)*. 

The final point to be raised regards the Russell 
model and rectification. The results obtained with WINK8 
Suggest that there is a large reflection effect. Hence, 
one may ask whether or not the Russell model, as used in 
“the methods of Tabachnik and Kitamura, is really valid. 

To answer this, it should be remembered that the result 
obtained with WINK8 involved holding the secondary tempera- 
ture constant, and letting that of the primary vary. As 
was mentioned earlier, this would influence the magnitude 
Of the reflection etE~ect. One should, therefore, treat the 
WINKS result with some caution. Also, the Russell model is 
a first approximation, and serves tO provide a preliminary 
set of elements. During the eclipses, reflection is being 
ignored. Furthermore, to make use of rectification would 
not be correct, because one cannot reliably extend a fitted 
Giinve beyond. the Gomain Gf the Flt we lne sr ectit ication 
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A +A, cos@ +A,cos‘6) are fitted to the out-of-eclipse 
variation, and the resulting Fourier series is applied to 
all data points, both inside and outside eclipses. From 

the point of view of obtaining a set of preliminary elements, 
the use Of the Russell model tsiqustiafied. 

In general, it is not possible to establish a set of 
criteria for choosing the best elements. Therefore, a 
choice of a "best" set elements should be based on the 
observer's experience, and on any other available informa- 


tion regarding the particular binary system. 
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CHAPTER 8 


W DELPHINI 


Sel introcuct ion 


The second Star to be analyzed with the various 
methods or Tight .cirve analysis Gs W Delphini. This star 
has served as a “standard example" in earlier work on light 
cunveranalysis (see Russell (1912b), Aitken (19357 p. 192), 
Tew CLIG2 7p. 597)5) TabachnLk and Shull berg (1967). and 
Tsesevich (1971) ). The light curve of W Delphiny has) been 
well-observed by Wendell (1909, 1914), and more recently by 
Walter (1970). The data to be used in this chapter are 
taken from both of these sources, in particular, the pri- 
“mary eclipse data is taken from Irwin (1962), who used W 
Delphini to demonstrate the Russell-Merrill method (Irwin's 
data is from Wendell's light curve). 

W Delphini has been classified as an "Algol-type" 
{i.e., detached) system by Russell (1912b), although the 
results of Walter (1970) indicate the presence of mass 
transter and a Significant reflection eGirect. “The orbital 
period of W Delphini is 4.8061 days. The primary eclipse 
is a very deep total eclipse (flat minimum), the magnitude 
change being 2™5695. “The secondary (annular) eclipse is 
veby shallow and almestinon-exustent, | liwin (1962), p'.599) 
has computed the change in magnitude to be 0".039. The 
Deimaryeecclipse of W Delphinw is shown ain seigime Loe ein 


table 6.1, the vallies of the geometric elements, as found 
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FIGURE 19, LIGHT CURVE-W DELPHINI, 
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DY sEWwirie (D9 627. D. 590 laresci vers melt eneulan ee noted that 
thesvalue of x= 0.5 used for the limb darkening by. Trwin 
was) fOL Cemonstration purposes only. If the spectral types 
OF LNG Stars are taken ante account, a more realistic value 
as x=0.6. The spectral types, as given in the catalogue 
OP Bacten et al (1982) are Ale tor the primary, and G5LvV 
hor hey secondary, “Furthermore,. the, eccentricity €= 0.20, 
Ehesmassy sUneCE LONI (m)-—10.0127,) and asin . = 1.94 Pall km. 
The orbital elements of Irwin, “as gGivén in table 8.1, will 


be used aS a cOmparison with results obtained by other 


methods. 


So208 Uuighnt Cunve Analysis 


The results obtained with the methods of Tabachnik, 
Kitamura and Kopal appear in table 8.2, and the WINK8 
results in table 8.3. The lvght curve of Wendell was 
analyzed using the methods of Tabachnik, Kitamura, and 
Kopal, while Walter's light curve was analyzed with WINK8. 
(The observations were weighted according to the square root 
Of the number Of Observations Made aula spacticular phase). 
Tabachnik's method was programmed into a TI-59 programmable 
Galculator. Once again, the program LCEIT2 computed the 
quantities necessary for an analysis by Kitamura’s method. 
An analysis using Kopal"s method was carried out using the 
Kalman filter (program KAL) and the total eclipse program 
TOTAL. The Kalman filter was used because both light curves 


had telatavely few points (Wendell: 27 points, Walter: 59 
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Tapes Gl) WW Debohini — lrwin 3s fh lementes 


Ladiusso: large: star = Be = 052474 (k= 0) 2612 5)) 
radius) of small star = Co OFS 46 

b=" 85..07 
Onto polarcgees tar.) s— 7 = 0.9164 
Pignhevoeasmall= star —= L. = 0.0836 


limb darkening x = 0.6 


Table 8.2 The Geometric Elements of W Delphini 


Tabachnik's method (TI-59) 


i = 88°.04 
radius =ot large star = r= Ons PAPAS) 
radius Oe small “star = Lon O76 
iamb darkening x = 106 
L, = 0.9164 L. = 0.08360 


Kitamura 's method (LCFT2) 


i =F80. 

Lee 0.24 ry = Oe aul 

= 2505556 

= 0. = 0.08360 

Ly 0.9164 L., O08 
Kopal's method (KAL, TOTAL) 
radius of eclipsed star = r= OF, Was: 
radius of eclipsing star = ry = 0.241 

i = 86° 

L, = 0.9164 L, = 0208360 


(Note: L, = light Oi larger Sta, L, = light GieSma ote tai 
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Table 8.3 WINK8 - W Delphini 


i= 35° 
PACING TOLeDiliiaty star = r= Oa Ay 0 0016 
Padus OL Secondary star = ry = Oie4 0. O22 
L, = 020213 
Ly = 0.07868 2) 2 Sse 
Astrophysical parameters 
star v ay 5 (equator) “T(polar) “ogg polytrope 
A Orl4505 7) 021480 0-250" 960050" K 9623.70 *K ='420 52.0 
By, (0243198 8022432502 250' ~ 46992 02K A796. 0K = 450 BO: 
Model parameters 
a = 0.14863 kK. 21 ot K = 1.64284 
Jo 310 = 0.02741 J worm = 9° 92723 Je o1o = 9° 05742 
mass ratio=0.5 quadrature magnitude = 0.025 
Stare Clatpeicicy e Uy=x Uy, Vv W a b c 
A OF 0997. 67 0.0 0.60 0.0%—450 =150 "0.1486 071463" 0.1479 
B 09509 0202057" 0. 60F 070-4: 0 120: 9052555" =O 22450 #027598 
Luminosities 
star apparent normalized total (4T) normalized (4mT) 
A 005538 Ono 2a 2 0.1299 10-- 0.86606 
B 0.00473 0.07868 0.2009x101* 02 13394 


ratio 0.08540 0.15465 
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points). In such a case, least-squares Fourier analysis 
would not produce meaningful results. This was substan- 
tiated when program EB.FS was used to compute the moments. 
several trial runs produced situations in which some of the 
moments were negative, had large errors (greater than the 
moments themselves), or whose magnitude did not decrease 
with increasing values of m. On the other hand, the Kalman 
filter produced positive moments, whose magnitude did 
decrease with increasing values of m (see table 8.2). 
Unfortunately, the use of the Kalman filter precludes any 
error analysis. Moreover, the elements produced with the 
aealie: computed by KAIewere much move reasonable. “Other 
sets of moments (from EB.FS) produced values of the primary 
star's radius greater than one, which is an unreasonable 
value since the light curve 1s mot One Of an ~~ over—contact ® 
(i.e., Similar to W Ursae Majoris) system. The analysis 
Using WINKS tdid not produce any unusual results, ‘other "than 
ease 1t found both eclipses “Ee be “partial, racher than 
total and annular. However, such a result is not unreason- 
anole, giveneulat Only 59s points Were ised, angmelaw ee 
eciupseisrardrazing (i.e., gust bavely eCOtal) One = Hence, 
asitaki change 2n the inclinatiron couldvalter tne type 16: 
eclipse crom total to Darulal an seve Cie oer Cults ee tiene 
is arstonireicant reflection eifect von the yecondary star 
(29702)°*K) but«a negligible ettect on the primary. ‘The 
stellar semi-axes showed that both stars were spherical, 


Wlth the Secondary being slaghtly more polar-tlattened chan 
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the primary (see table 8.3). The complete analysis required 


bive .terations. 


Geos conc luszons 


In concluding the analysis of W-Déelphini, there are 
some points worth discussing. As with HS Herculis, there 
are differences between the different, sets.of «elements. 
However, the difterences are between the methods uSing the 
Russel!) model and those mot ising this models. Theseléments 
obtained with Tabachnik's and Kitamura's methods produced 
nearnly.identical results (except for the radius of the 
smaller star), while Kopal's method and WINK8 produced 
identical sets of elements. Once again, one is faced with 
the choice of choosing the best set of elements, but given 
the more realistic models employed by Kopal's method and 
WINK8, one would tend to believe the sets of elements 
obtained with these methods. It should be noted that a 
ALEreLrentilal Conrectilons, pregram was not wsed, Simee) the 
system is well-separated, and also for the purposes of a 
direct comparison between the methods. One can also argue 
that the elements. obtained via Kopal's method. are closest 
to those found by WINK8. Any further differences between 
the various sets of elements should be attributed to the 
computational methods employed. 

The analysis just presented reveals an important 
pointe negardaing Kopal’s method. anethe previous section, 


4+ was mentioned that the method of least-squares Fourier 
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analysis produced spurious results, making the use of the 
Kalman filter a necessity. This Situation precluded any 
SOnG Of @rror analysis. ‘Hence; Wwe one wishes to obtain a 
set of elements with error estimates using Kopal's method, 
one must have a large number of data points available. 
THLS Situation also demonstrates the Utility of the Kadman 
filter and its ability to produce good least-squares esti- 
mates of the moments using relatively few (27) points. 
Therefore, the Kalman filter is well-suited to situations 
in which the light curve consists of relatively few data 
PCaintS. = ne Oppesice would) be yerue: ©Oreteast-—squares 
Fourier analysis. 

As a final note, this analysis of W Delphini shows 
that many of the eclipsing binaries analyzed using Russell- 
type methods would bear a second look, both from the point 
of view of light curve analysis, and from an observational 
point of view. Systems such as W Delphini should be re- 
analyzed using modern methods (Kopal or WINK8), which would 
use well-observed (i.e., several hundred data points) light 
curves as input data. One could then compute reliable sets 
of elements, with good error estimates. In this regard, 
this analysis of W Delphini is somewhat original Since this 
System has not (to the author's Knowledge) been analyzed 
with any method more complicated than Tabachnik's method. 
Rather than reatftirming old sesulius,, anv analysis such as 
the sone just presented provides greater insight into the 


phenomenon of eclipsing binaries. 
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CHAPTER 9 
HD 219634 


eM Mek sielexe Uhvey aalove 

ine wlast eclipsing binary tO Dervconsidered ianechis 
survey of light curve analysis methods is the recently dis- 
covered system of HD 219634. The system has been observed 
photometrically by Gulliver, Hube and Lowe (1982), but has 
been monitored only intermittently since the discovery of 
its eclipsing nature. Hube has also made several spectro- 
graphic DASE INS. with which a spectroscopic solution 
has been computed (not yet published). The spectroscopic 
solution revealed a mass function of 0.16 solar masses and 
a period Of 233912 days. The large value of the mass func 
tion suggested that HD 219634 might be an eclipsing binary. 
Gulliver et al have deduced a spectral type of BOVn for 
HD 219634. Since the secondary was not detected spectros- 
copically, this classification would refer to the primary 
Stare) Lt should also be noted that BHD 219634" tsa) posstbie 
Optical counterpart to the X-ray source 4U23416 4 61 (see 
Horman tet al (1978) )~ 

The partial lightmcunye sO beuuSsed sll lilcMaial yous 
Was Obtained by- GubLiverwate Miele Peakwi nes Od eus tinge cist e— 
bentraleUBVY photometry. | The AV dightvcunye vappears an 
Brguresl Oe Theslight curve Naseseveralwnoelceab leu toacurese 
Theweclipses appear tO be partial, since the minima yare 


Earniy sharp. | Secondly, senere Usa) noticeable lack ot 
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points outside the eclipses, a feature which automatically 
rules out. least-squares Fourier, analysis, in. Kopad 4s, method. 
Hence, the, Kalman, filter must. be msed.. The oniginalsslaght 
curve of Gulliver et al was asymmetric in the sense that the 
shoulders Obwtheslight curve, (points of maxterna ls tangency) 
were not, level »particularly at primary minimum... The asym— 
metry amounted to 01402 in the AV wight icurve.- The: Light 
curve shown in figure 19 has. this variation removed (all 
points with 0.0 < 6 < 0.0930 were shifted up by 0.02). 
Cuddiver et al, have attributed. thisevaria tion wo.2 mossibly 
Den cCOnstant, Compara son. stan (HD 220057) or {toian intrinsic 
Variation. in HD. 219634 itself... lt. should also iwbe noted: Ehat 
the light variation outside the eclipses is not constant, 

SO One would expect the components to be somewhat oblate. 
The observations to be used in the light curve analysis 


appear in appendix 5. 


Ok ohie CuLnvey Amat Sis 


The two previous chapters have demonstrated the rela- 
tive performance of the various methods of light curve 
analysis. HD 219634 ,allows these. methods to be applied co 
a new, and previously unstudied system. 

Since: the. light curve: of HDviloo24 shows ebreccs: due 
to oblateness and reflection, the Russell model is techni- 
cally not applicable. However, the method of Kitamura was 
used since the Fourier transform operation minimizes irre- 


Gilearittes: inethe light curve, andes; herefore preferable 
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to Tabachnik's method. Kopal's method was used for the 

same reason. WINK8, which models both reflection and oblate- 
ness in a realistic manner, gavé the best results. The 
results from the methods of Kopal and Kitamura appear in 
Saplewo == lane WwiINK@e results appear in table orc. sed 
methods produced an inclination close to 69° degrees. How- 
ever, the agreement between the various radii is not good. 
In all cases, the radius of the smaller star lies between 
0.1 and 0.2, while the larger radius lies between 0.12 and 
0.36. A similar situation exists for the relative lumino- 
Sities. Since WINK8 uses a more realistic model, the 
results so obtained are probably closer to the truth. The 
value of k=0.799 obtained by WINK8 implies that the primary 
Gcilipsewts altransit (ise. , annular). 9)in using Kkopal’s 
method, it was found that the primary eclipse was an occul- 
tation. This was determined through the use of the Russell- 
Merrill y-functions (see Appendix 3). No such assumption 
was required with Kitamura's method. The values of the 
eharacteristic Lunction '» indicated that the primary eclapse 
isvan occultation, agreeing with" Kopal’ s method.” Further 
analysis to resolve this disagreement between the Fourier 
transform methods and WINK8 is most certainly required. 

Some other interesting results obtained with WINK8 were that 
Che warpgermestar 1s also the brignitem one, and that, Unere yrs 
a ’slagnificant reflection effect on che secondary” (Al 9250 "K). 
Both stars are also significantly oblate. The temperatures 


Obtarned) by WINKS "did not agree with the spectral classifi— 
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Gabon sot (Gulliver et aly;iwhich was BOVn )(29000°) sfor the 
Primary fand Be (12000°K) for ehestsecondary. this arises 
from the fact that WINK8 was allowed to use its default 
parameters in starting the solution. The default values of 
the. temperatures are both 0000°K. s Thus problem wath the 
temperatures will be discussed further in the conclusions 
of this chapter. The discrepancies in the radii and rela- 


tive luminosities displayed in the methods of Kopal and 


Kitamura can be attributed to the computational method used. 


in the case of Kopal's method, an initial estimate cf k was 
Eequired, sand the resulting elements reflect the mather 
large uncertainty in the value of k. Unfortunately, the 
method for determining an improved value of k led to wildly 
erroneous values, so only an initial approximation to the 
elements was possible. In the case of Kitamura's method, 
Pe) Was, possibile «bO, Match «the: characteristic GunctiLons only 
for the case r= 0.36 and r= 0.13. Thus the problem lies 


Witch “any wncertainties: that) the characteristic munctions 


might have. However, the radii so obtained are closer to 


those found by WINK8 than are the radii from Kopal's method. 


OnGe. again, similar! remarks would apply to the) relative 
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Ore nCOnclLusions 
ineconucluding \enis analysis ©: ED 71 [ola eenere, ake 
some points to be discussed and some observations to be 


mades ini the previous: Sections 1c was mentioned that the 
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temperatures found by WINK8 did not agree with spectral 
classification of HD 219634. A possible remedy would be to 
fix both temperatures at the values determined from the 
spectral "classtilcationm and co computermanotier isolucion:. 
The temperatures would affect the solution through the limb 
darkening (which was held fixed at 0.6, according to the 
stars' spectral types) and through the value of the central 
intensity, I(0), Since it is only these parameters which are 
directly affected by the temperature. However, one should 
not rule out more subtle connections between the tempera- 
tures of the stars and other parameters. Since a change in 
limb darkening can mimic a change in the radius, the limb 
darkening was held constant. Hence, the only direct effect 
that a change in temperature would have would be on I1(0), 
which ultimately determines the luminosity of the star. 
This would obviously have an effect on the radius of the 
star, although the magnitude of such an effect could not be 
easily estimated. Hence, a further avenue of research would 
be to investigate the change in the radius induced by a 
change in the temperature. Clearly, the results obtained 
with Kopal's method require some comment. The problems lie 
in the computational method used. This method is outlined 
in Section 4.9 of Chapter 4. Using this approach, it seems 
EOmDeECULECeCdLeticuliey 1 fanots IMpOssible, Sco proceed) bOoa 
more refined solution. The problem with the present method 
Gf solution in the case of partial eclipses is that an 


initial value of kK must be. speciiied., A, better approach 
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would mean abandoning the method of Section 4.9 and to use 
a more fundamental approach. Kopal outlines such a method 
imnis 1979 book (pp. 175-179) Wathout going anto: creat 

detail, the method involves finding the roots of simultan- 


eous non-linear equations 


(aigad Uelg ores) hv and g,(a,c_) = 
2 Oo Asa, 4 fe) AoA 
ah a si 
where By ae and eo = —— : 
aba, a 1 


The right hand sides of the first set of equations are de- 
termined Geom the observations. The left hand sides are 
then computed using the series expansion, equation (4.48) 
of Chapter 4. The problem is then an iterative one for 
determining a and Co: 

Finally, the question arises regarding the best set 
of elements. It was mentioned in the previous section that 
the set of elements that are closest to the truth are those 
found by WINK8. Since WINK8 fits the whole light curve, 
and accounts for both oblateness and reflection, one would 
expect the WINKS elements to be the best. The only other 
solution that comes anywhere close is the set of elements 
from Kitamura's method. With regard to WINK8, the problem 
of the discrepant temperatures should not be forgotten. 
Clearly, HD)219634 mMer1ts further investigation, boOel Obser— 
Vacionally and analytreally.. Ene wnums>er son Tout —oE—ee lise 
observations is too small to reliably model the second- 


Order effects, and the source(s) of the displacements of 
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Ehe two branches, of each eclipse 2s still noe known. in 
Spite of their small number and (possibly) large uncertain- 
ties, it is comforting to note that some sense can be made 


of the observations of this system. 
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Table 9.1 HD 219634 Orbital Elements 
Kitamura 
I= 160" 
re (small star), = ,0f12 


G Charge tstar jy = 0.38 


b 


ith (pramary) ==" 0208106 


Le (secondary) = 0.9189 


x = 0.6 

Kopal 
i = 0-63. 4.4 
ry (eclipsed star) = 0.168 
ro (eclipsing start) = 0-122 
Ly (eclipsed star) = 0.2920 
L. (eclipsing star) = 0.7080 
x= 046 
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Table 9.2 HD 219634 Orbital Elements from WINKS 


eet ishe Bue {Os Wanel WIS ashe 


magnitude at quadrature = -0.429+0.001 
Ty = 10000°K (fixed) 

T, eo A Osher 2. lke 

r, (large Star) = 0249 2 0.004 


r., (small Siar)? = SOR or 250 0s 
Ke sO O98 O04 2 
x= 06 (fixed fom both Stars) 
)(0-c)* = 0.00793 


Astrophysical Parameters 


poly- 
star Vv aa 8 Tlequator,. 7h) “P(polary,,°h)s) log ig) 2repe 
A Q25239. 0225245 50-25 10000.0 LOIS5, 81 4.0 cal 
B G2NS827 O21 883.5 0225 902726 JisIeLs 4.0 5a 
Model Parameters 
a= 0.26054 ig 0. 73230 So Oe 4:5 9;7 Je5Q9 7 0.78784 
ioe = 0.78873 Dea = 0.670 50 ome ete 
star ellipticity C u, (=x) us Vv W a b @ 
A 0.9766 0.0076 (026> 020 (—450= 1207 90-2605" 9022545" 022502 
B 0.9901 0.0033 036° O20" —-470 (120 1021908" Or1ess Soe 77 
Luminosities 
Siac ia apparent normalized otal, C47) normalized (47) 
A 0.16460 0.69642 0.4482 x No Onwi26 97 
B 0.07180 0. S0S5c OkebGs 3x hee 0227563 
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CHAPTER 10 


CONCLUSIONS AND COMMENTS 


After having applied the four methods of light curve 
analysis to three different eclipsing binaries, a number of 
questions quite naturally arise. To begin with, how do the 
various methods of light curve analysis compare? How physi- 
cally realistic are the various models underlying these 
methods? 

The analyses in the previous three chapters have 
clearly established that Wood's model is the "best" since 
the elements obtained are the result of a fit to the whole 
light curve, rather than to one half of an eclipse. Agree- 
ment with published results confirms that WINK8 is operating 
correctly. However, Kopal's frequency domain method per- 
formed just as well as WINK8, despite the fact that a spher- 
ical model was used. This is particularly evident in the 
analysis of W Delphini. One can therefore say that Kopal's 
method can produce results as good as those found using 
WINK8, but with relatively little computation. However, 
Kopal's method does not provide the detailed description 
that Wood's model does. For instance, one cannot estimate 
the magnitude of the rerlection Gliect or Compute the separ— 
ate Ssemiaxes a, Db, and C in KOpalys mothod = Also, Kopal™s 
method does not really allow one to compute a theoretical 
lightecusye, abart: .LOMm COMPUtING s0eecxDilcivLy, samprocess 


Which 1s quite tedious 22 the stars on] particular, system 
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aver signibicantly distorted ys.) Wath tehese Vimitc, however, 
one can regard the methods of Wood and Kopal as being equi- 
valent. The two methods incorporating the Russell model, 
namely the methods of Tabachnik and Kitamura, show the 
@ppBRoximate mature of the modelys anduhence tts: dnabiiity to 
represent systems which do not have spherical components, 
or those in which reflection may be important. Of these 
two methods, the more desirable of the two is that of 
Kitamura, since the Fourier transform operation acts as a 
smoothing process. The success of Kitamura's method in this 
regard ¥s,evident in the case of HD 219634, even though the 
elements were slightly different from those found by WINK8. 
A disadvantage that both of the Russell-type methods have 
is that they are indirect. One must compute various auxi- 
liary functions and then either compare these with theore- 
tical computations or use «them inva further computation, 
which finally produces the set of elements. On the other 
hand, WINK8 and Kopal's method proceed to the solution more 
directly. In the case of WINK8, the residuals AI(0-C) have 
a direct bearing on the resultant elements. In Kopal's 
method, the moments of the light curve enter the equations 
for the elements directly. 

As far as realism is concerned, the points made in the 
Previous paragraph apply= The models of Wood and» Kopal® are 
Cleaniy ssuperion to: therkussel) Iamedei@scincert ies physical 
Properties of the Stars enter Enevanalycis more: directly. 


The Russell model is really little more than a geometric 
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one. With the models of Wood and Kopal, one can represent 
stellar oblateness and reflection rigorously, © rather than 
resorting to the questionable approach of rectification. 
Also, the models of Wood and Kopal benefit from our current 
knowledge of stellar atmospheres, especially with regard to 
limb darkening. With regard to the reflection effect in 
particular, the models of Kopal and Wood handle it in a 
Satisfactory way. WINK8 models the effect directly, whereas 
Kopal's method allows one to filter out the reflection 
effect. 

The. analyses of the previous chapters are but prelim- 
inary Ones, “particularly: in’ ther cases of UD 219634. "Arter 
such initial analyses, one can go on to refine the solution, 
investigating any interrelationships that may exist between 
the various model parameters. A case in point is the "dis- 
crepant temperatures" problem inHD 219634. Clearly, further 
work on this system is were ceee If the elements of a 
particular eclipsing binary are reasonably well determined, 
some other avenues of research are to determine the eccen- 
tricity of the relative orbit; the absolute dimensions of 
the system (using spectroscopic data), and to investigate 
the spOSSTDibLaty Of apsidalumotion. “Theslast 16m is parti 
Cularly important with regard to relativistic effects and 
stellar evolution. Hence, one can gauge the importance of 
mass transfer. An independent approach would be to examine 
the spectrum Of a system,  whicm provides ene with some idea 


of the rotational velocities of the two stars. There are 
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also Other ertects that can be investigated in the spectra 
Of eclipsing binaries. The analyses of HS Herculis and W 
Delphini show that these systems are worth a second look, 
Since the newer approaches to light curve analysis provide 
Greater Insight, into Ehe nature of veclipsing binaries, 
rather than reaffirming older results. Many eclipsing 
binaries analyzed with the Russell model would certainly 
benefit from more detailed consideration using one of the 
modern techniques. 

In conclusion then, it can safely be said that modern 
methods of light curve analysis provide greater insight 
into the nature of eclipsing binary stars, since the models 
employed and the methods of analysis used are based on more 
realistic physical models. However, a "best set" of ele- 
ments is not easy to obtain since light curve analysis 
methods are not "black boxes" which produce results indepen- 
dent of the user's judgement. Light curve analysis by its 
very nature is not a simple problem, so an "experimental" 
approach is required to really come to an understanding of 


the physics involved in eclipsing binary Stars. 
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Appendix 1 Russell Model Programs and Rectification 


Al.1 Computer Programs 

There are several programs which use Russell's model 
in the analysis of eclipse light curves. The FORTRAN pro- 
grams LINE and LINE2 use Tabachnik's method. LINE takes a 
Value of k and the values of @ and p aS input data, After 
choosing an appropriate value of k, one must use the Russell- 
Merrill tables (see Russell and Merrill (1952)) to find the 
values of p corresponding to the values of a computed from 
the light curve. LINE does all other computations. LINE2 
is ae LINE, but takes as anput values of, p corres— 
ponding to two values of k (maximum difference of 0.1 ink 
allowed). LINE2 then scans all k values between the chosen 
lim2ts, “and prints out the corresponding values of yee e 
The value of k for which dna) re is a minimum should be 
taken as the solution. Both LINE and LINE2 make use of an 
imeecer array known ‘as INDIC. If, ENDIG(1)—1, then a par 
ticular observation is included ain the least-squares £it. 
An Observation 2s atte we INDIC (CE) 4 ib Wwalues for INDIG 
follow those of 6 and p; 1s2e5, Cie Pas L, tor all observa= 
tions). The programs are initiated by the commands 

Si ea 

¢$R *FORTG SCARDS = LINE T=1.0 

(Computer will respond with a message.) 

$R -LOAD# 5 = (A DATA ELL GE) 6 =-P 


(another message.) 
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Torsee che results; lise the output 21 le" —P. 2 Al FPORTRAN 
PLogrvamsmare weed “in this fashion. A program using 
DTabachnik “s' method is also Given for the’ TI—59) calculator. 
The instructions are self-explanatory. The programs ZAPP.DC 
and” DLGHT "ace -ditfierential correction programs. ~ ZAPP. DC is 
a WATFIV program that uses values of the derivatives from 
Irwin's tables (see Irwin (1947)). The generalized linear 
least-squares program used in ZAPP.DC has been adapted from 
the program published by Malin et al (1982). The values 

of the derivatives are found by interpolating in Irwin's 
tables. One need only enter the value of a, the values of 
the derivatives, the value of A&(0-C), the weight (usually 
1.0), and the phase, for each observation. ZAPP.DC then 
formulates the equations of condition, and solves them for 
9! and WSR) The value 


Of the= lamb darkening, which is required for the solution; 


the corrections Es oye Ari: ‘Nie 


is’ noe corrected ~ The program isting ancludes a “sample 
input. The standard deviations of the corrections are also 
PEeinced out, as is the value of aise Avcample output 


follows the program listing. ZAPP.DC is initiated with the 
command 

$R *WATFIV SCARDS = ZAPP.DC 
Al sresuites are prianted out on Ehe "cerminea y. “Program 
LIGHT uses an entirely different approach. All derivatives 
ace computed numerically,. and therrequrred ¢—-tunce tions are 
computed using the polynomial approximations of Fliegel and 


Wilson (1968). The least-squares procedure used is known 
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as the Levenberg-Marquardt algorithm (see, for example, 
Gall, Mizray and Wright (lO Gi opel 36-137). SAP pRogram, 
Known as CURFIT, which incorporates this algorithm has 
been published by Bevington (1969, pp. 237-239). This pro- 
gram, aS well as other required subroutines have been 
incorporated into LIGHT, UGIGHT requires only the input 
light curve (eclipse part only), as well as starting values 
for the elements and step lengths to be used in the solu- 
tion. Due to its large size, LIGHT is not reproduced here, 
butea sample? run? and data fillet are given. "LIGHT i1s*>initia— 
ted in the same way as LINE and LINE2, but one must use the 
following sequence to run the program: 

$R - LOAD# 4=COEF 5= (A DATA FILE) G= —P 
File COEF contains coefficients for the polynomial appro- 
ximation to @. It should be noted that LIGHT cannot handle 
the annular phases of annular eclipses because an approxi- 
mation to a for such phases does not yet exist. Such an 
approximation will be incorporated in future versions of 
bt Gh. 

There are a few subtleties to be aware of when 
using L.tGHl. VOne should not use = large seep lenges 
("02.4) Since this can produce: wildly erroneous results: 
Hence, “allustep. tengths should besne Wargerns than. 0701, a 
Goode typical value bende 0.001 gene Contr dai e-OnceG nae 
entwval corrector sucheas, ZAPP<DG,) LIGHT prints ouc tne new 
set of elements followed by their standard deviations. 
Corrections, are Nol Printed, OL, ) DULZa 1S. Otelnpu en para— 
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Al. 2 An Example of Rectirtication 

ihre method of rectification used 1s the Graphical 
metnod Of Russell “and Merrill. (1952, 0g..54)>, and is) des-— 
eceibed) completely Cherein.. (in, order to rectirty 4 light 


Gurve, one fits the Fourier series 


L(G) =AL+A cos 9 + A.cos 26 + B.siné@ + B.sin 26 


1 Z ui 2 


to the light variation outside eclipses. It should be 
noted that the presence of the last two terms cannot be 
justified physically. One is also assuming higher terms 
to be negligible, which is not usually the case. The 
method itself is quite simple: one reads values of &(8) 
from the light curve at some convenient interval in 6 (say 
LOs Opis. scad), letting a, b,c, dad mepresent: 2.(6) mat) oy 
Le0r— 0 7) 60 oy and.) 300 .—0.. By COMbinang a,b, Cpnand Gio 


various ways, one finds that 


—~(at+tb+ctd) = A + A, cos 26 
Le Si = Bere = A. cos 6 
4 dl 
Di Gown Be sine: 
4 A 
Dees oe ey Ba siniae 
4 2 


Mheretore, by plotting the Lette hand side of these equa-— 
ions against cos 20, cos. oi, Sin @, ana Sin 20> One can 
find the coefficients from the slopes (and intercept in 
Ener tare. equation cnly om the, lines Any pointes saappen= 
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so one may effectively exclude in-eclipse observations, 
Using this criterion. The reetiiication, proceaqure is then 
Cabried Out using! equattonse2. 25) and Ore) On chapter... 
To illustrate the procedure, observations of HS 
Herculis outside eclipse will be analyzed. To facilitate 
the computation, only a small number of observations will 
be used. In general, it would be advisable to use all 
available observations. The procedure described above can 
be performed with relative ease using a programmable calcu- 


,A B, and B. may then be 


L Doe eel! 2 


found using a linear least-squares routine. The results 


PAtoOr. The coefficients Ass A 


of the calculation are presented in table Al.1, along with 


those Of Hall and Hubbard (4971) 


Table ALi: 
present work Hall and Hubbard 

AG OF 9376 Sr Sreh0) 
Ay -0.01645 =0.015 
A. a= (nhs Seon oO 00s 

. Or0 
By Or0 
B O58 OmA0 
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COMPUTATION OF ECLIPSING BINARY ORBITAL 148 


1 e 

2 C ELEMENTS BY TABACHNIK*S METHOD. 

3 C PROGRAM ‘LINE*e BY DeHOLMGRENs JANe8&' 34.6 
4 Cc 

5 IMPLICIT REAL (A—H.0-Z) 

6 REAL INCL 

7% INTEGER INDIC(20) 

a DIMENSION T(50)+P(50)+sX(50) s¥(50) 
9 © READ IN DATA (THETA.P)- 

10 READ(5S+1)NOBS+RK 

14 1 FORMAT(135F14.7) 

7 READCS45 22.) (1 CI) se U1 es tND1e@l).) 1=15 NOES) 
13 2 FORMAT (2F14.67,.13) 

14 C FORM X AND Ye 

15 DO 3 I=1,NO8S 

16 MOTI =SINCT C1) eSIENGHICIOD 

17 BEV CN=C1 6 tek te Ci) x Clete ReseCr)) 
18 C DO LEAST=SQUARES FIT 

19 SY=06 

20 SX=06 

21 SY2=06 

22 SX2=0-6 

23 SXY=06 

24 BOP4)J=1, Noss 
25 meet (IND TCO s)6Ne. io) eGo toma 

26 »  SxX=SX#X(J) 
2 SY=SY#Y(J) 

28 SX2=SX2+X( J) *X(J) 

29 SY 2=SY2+Y(J)*Y(J) 

30 SXY=SXY#X( J) *Y( J) 
31 4 CONTINUE 

32 C DETERMINE PARAMETERS. 

33 AA=SXY—( SX*SY/NOBS) . 
34 3R=SxX2-SY2+( (SY*¥SY-SX*SX )/NOBS ) 
35 DSC=BB *BB/ (40 *AAXAA) +1 

36 Te COSC fle em 0 NDSC=—0SC 

37 A=-BB/ (26*AA)+SQRT(DSC) 

38 B=USY=A*SX)/NOBS 

39 VE CAslo hs. OL) A=—A 

40 Te CS) 6121) 10.) S=—8 

41 TI=SQRT(A/B) 

42 INCL=ATAN(TI) 

43 R1=1e/SQRT(A+t8) 

44 R2=RK*¥R1 

45 WRITE (655) Ri shes 1 NCE 
46 Sy CORMATN et RI= "5h 14s 7s Xe Ree 'sF 147s eV INGCL=".hl 467) 
47 © FINDIESUM SQ. RESID<. 1. 

48 S2=06 

49 DO 6 K=1,NOBS 

50 6 S2=S2+ (A¥X(K)-Y(K) +B) * (AKX(K)-Y(K) 4B) 
51 S2=S2/(166t+A*A) 

52 WRITE (657) S2 

53 TFOEORMATCOSUM SQe RESID. ="*sF14e7) 
54 stop 

55 END 
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COMPUTATIGN OF ECLIPSING 
ELEMENTS Gy TABACHNIK*'S: METHGCDs 
PROGRAM ®LINE2%*.« BY DeHGLMGRENs JANe22' 846 


IMPCICIT REAL (A—H,O—Z) 
Rev INK 
PNT EGER INOTC C29) 


DIMENS TON 1450) 2P(5S0) 5X1503467Y(5S0)sPLO( SCG) .,PHT(50) 


C teeesyS) UN) IEA ON UCUinre yr 5 = pA 
READ(541)NOBSe NKeRKLUs KKHI ,DEL 


1 


2 


FORMAT (23353F1427) 


eS OND) (CS) 5/2 7) CT Gin) ol OXelis > Ht Cis Nw Ce Gi), 


FURMAT (3F1427,13) 


© Fer A AND® Ye 


3 


4 


RK=RKLCO 
DO 100 L=1,NK 


CALL INTER( PLO, PHI sPsRKsRKLOGsRKHI »NORS) 


DO 3 I=1,NO8S 
XCTI=SIN(TCI)I*SINOTI(T) 


EiNARY CREITAL 


) 


YCLTI=HCLetRK*P(1))*( 12 +RKXP( 1)? 
CY DOVE EASA—SOUARES nia 


SY=0e 

SS 85 

SY2=0. 

SX2=06 
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DO 4 J=1s.NO8S 
SCIEN D) RES Ce) GING sl) IGG) Ue 
SX=SX+X{ J) 

SH\(S Siar if CS) 2) 
SX2=SX2+K(J)*X( J) 
SY¥2= Sve y (9) = Y (32 
SX Y= Sn YR ee Ke) 
CONTINUE 


C DETERMINE PARAMETERS» 
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AA=SKY-(SX*SY/NOBS ) 
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BB=SX2—SY2+(USY*SY-SX*SX)/NGBS) 


OSC=BB*BB/(4 e*AAXAA) +16 
Te WOSEG yells) Oe USS ——DSG 
A=—-B8B/(2+¢*AA)+SQRT(DOSC) 
B={SY-A*SX)/NOBS 

TF (AeL Te00) A=—-A 

IPC elie 06) B—=—5 
TI=SQRT(A/B) 
INCL=ATAN(TI) 
RL=1e/SQRT(AtS) 
R2=RK¥RI1 
WRITECGs5IR1 sR 251 NCL 


t= 1 SMOKES 
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WRITE(6412) RK 
FORMAT (1X5 "*K=* »FI4e7) 


C FIND SUM SQe RESID. -« 


100 


Sic=Ole 
DO 6 K=1+sNO8S 


S2= 52+ (AX K)—Y¥(K 748) FCAEX( KD -—Y¥(K) 48) 


S$2=S2/(1¢etA*A) 
WRITE(6s7)S2 
FORMAT(* SUM SQe RESID. 
RK=RK+DEL 
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ol END 

62 © LINEAR INTERPOLATION, SUBROUTINE. 

oS SUBROUTINE INTER(PLOsPHI sP 5 RK»e RKLU ok KHI ps NUS ) 

64 HEE Me tty MeiNk tase OS <4 y 

Ey5) DIMENSION PLO(SG)5 Pha (50) sPi5G) 

66 DG 1 IT=1,NU68S 

67 1 PCIJ=PLO(41)4(P RIT 1) -PLO(1 ) 2 F1RR-RLOD/ URKHIM-RELO) 
68 RETURN 

69 END 
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PARTITIONING (OP 17) A Z 9 5 2) LIBRARY MODULE a PRINTER rs CARDS il 
PARTITION (Op 17) MODULE ENFICHABL eG « =—————— I MPRIMANT ES = —— = = CARTES 


PROGRAM DESCRIPTION ® DESCRIPTION DU PROGRAMME 


Program solves for the elements of an eclipsing binary (fx 
in thevcase of complete eclipses (te., total or annular). 


grt) 


Method used is a variation of Kopal's method due to V.M. Tabachnik, 
which fits a line of the form y=ax+b, where y=(l+kp)4 x=sin’ , 
ae 4/7 and b=cos7i/r7 (kK=rg/r,)). p is the geometrical depth, 


and is used as input data (use Russell-Merrill tables of a(k,p)). 


asf 
Elements are given by: tan i=a/b,r= (a+b) #x=kr, (k used as input) 


USER INSTRUCTION ® MODE D'EMPLOI 


STEP ENTER PRESS DISPLAY 
PROCEDURE INTRODUIRE APPUYER SUR AFFICHAGE 


Load program 


inieiali ze. E 0% 
Input k k r/s intermediate 
result 

EADitedata (05,0) @ where sO ii sain 
radians and p) is the corresponding 7 = ae a 
geometrical depth. D, we ra ee 

J points so far 
Calculate elemnts ("GO")- B ‘ 
first element computed is i. i (degrees) 
Display x r/s ry 
Display xc, r/s oe 
Fxamine correlation coefficient (CC) VE cc 
as a check on the value of k used 
(tome k close) to the correct ‘value, 
thesec. should: be. close, to 2) 
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TI PROGRAMMABLE 153 


TITLE On 
TIRESE eRe Se Mode |e ee ne pace Ok 4 
PROGRAMMER D. Holmgren May 19,1983 PROGRAM RECORD 

PROG RA EUR errr ee eee a) A eee es eee FICHE PROGRAMME 

PARTITIONING (OP 17) 

partition (or 17) L472 19,58 | SCC eta ee ee ee ee ee ee eae 


PROGRAM DESCRIPTION ® DESCRIPTION DU PROGRAMME 


Example:W 'Delphini - total eclipse. k is estimated from: 
k=O <6") /(6"46" = 10,5779 
where 9'=phase at first contact, @"=phase at second contact. See 


sheet of observations which follows program description. Results: 


k=0.5779: i=84.45(degrees) ca) 526s THe 3.40 
G=0.248 r=0.256 
y 0-144 r=0.135 (published) 


USER INSTRUCTION © MODE D’EMPLOI 


STEP ENTER PRESS DISPLAY 
SEQUENCE PROCEDURE INTRODUIRE APPUYER SUR AFFICHAGE 


Load program al il 
Initialize 
Enter k Onow a9 rsp IAS) 
Enter data 9=0.394 O. 
p=1.00 ies 
6=0.3304 A 1.1473 i 
etc., until all observations p=0.64 r/s De 
have been entered. ‘J 
Calculacre a B 84.45 
Calculate r, Es 0.248 
Calculate Ly ys 0.144 
Check Cé r/s 0299974. 22. 
Solution complete 
Ref. EClipsing Variable Stars, 
Chas, Edited by V.P. Tsesevich, 
John Wiley & Sons, 1973. 
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PROGRAM RECORD 
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PARTITIONING (OP 17) 4 7 9 5 9 LIBRARY MODULE os PRINTER CARDS q. 
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PROGRAM DESCRIPTION ® DESCRIPTION DU PROGRAMME 


Solves the depth equation using an iterative procedure. Input 


data is the timp darkening x (larger :ctar), Tay, A=00) Sak internal 


tangency of annular eclipse), A;depth of total eclipse (or value of £40) 


at internal tangency of total eclipse), and an initial estimate of k. 


The equations are: ¥ (k,-1)=(3 (1-xp) +2xp (11 (k) /k2) ) / (3-xp) 
eS Spy. 1t(k)=2_ (3sin-1/k-(3-4k) (1+2k) /k(I=k) ) 


Si 
USER INSTRUCTION ® MODE D'EMPLOI 


STEP ENTER PRESS DISPLAY 
PROCEDURE INTRODUIRE APPUYER SUR AFFICHAGE 


ha Y(k 


Load Program 
il Enter input data (xp,l-\4,4,,k) Bde SO 700 Xp 
(Anhitial estimate of 1) i 5 sta Ol t= 
ra STO 02 a 
k STOO k 
Zz Set radian mode - 2ND RAD 
3) Begin iteration r/s C flashing 
4 See new approximate k new k 
5 Continue iteration (as many times | 
as necessary to obtain convergence) cs 
*Requires only ~four iterations r/s finale 


with a good initial estimate. 


A bad estimate will extend the 
PLOCESSE bye 2 OL 3S aberrations. 
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PROGRAM RECORD 

PROGRAMMER 

PROGRAMMEUR Eo peO mg pene Slee tars DATE Man 27784 FICHE PROGRAMME 

tue veg) | Cah SE TET cree et ee a ee 


PROGRAM DESCRIPTION @ DESCRIPTION DU PROGRAMME 


Example: HS Herculis 
Xp=0.6, Ap=0.705, so 1=-’Ap,=0.295, Na-O0.9) Alnpira ly k=0n5 


final, k=0. 549410 (4 aterations) 


USER INSTRUCTION © MODE D’EMPLOI 


STEP ENTER PRESS DISPLAY 
PROCEDURE INTRODUIRE APPUYER SUR AFFICHAGE 


Load Deore 


i Enter data 20 OS 6S70 1100 0.6 
1l-dApH= 
022950. 295870 O01 OF 295 
Ne=0.9 (029 (ISTO102 O28) 
Ki ngeeoe> Uno SLO |03 es 
2 Set to radian mode 2nd |RAD 
3 Begin iteration r/s 1 
4 New k Oro 2.08 
5 Continue iteration r/s 
Next k value Os SIRES). 
GECc . 
Fourth iteration 0.549410. 
(No change in subsequent 
iterations) 
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PROGRAM *ZAPP.DC* FOR THE LEAST-SQUARES SOLUTIGN OF 

M LINEAR EQUATIONS IN N UNKNGWNS USING NGRMAL EQUATIGNS. 

THIS PROGRAM FOR DIFFERENTIAL COKRECTIONS USING IRWIN'S TABLES. 
BY Oe HOLMGRENs Oli el4*83e “SLIGHTLY MODIFIED — OCT. 20°33 


REAL TlseT2eWT ePV »SDsR1i aR2-5 INCL s RK es AL sP sDARK »sTAUsXTAU ePI 
INTEGER NTYPE 

DIMENSTO NBA CSAs 5 GSS cls e Ca 

DATA B43321¥* 040%, I1IP/B81¥*1/S 

PI=3.14159 


SPECIFY THE NUMBER OF EQUATIGNS (M) AND THE NUMBER OF 
UNKNOWNS (N).- 


READs MsNoNTYPE 
READsR1sR2SeRK»y INCL sAL »DARK 
INCL=INCL*¥PI/1 80-2 


etn Slaps IMs) 26i6) ee) he 
NP=N+#1 
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Appendix 2 Program for Kitamura's Method 


The program LCFT2 uses the procedure outlined in the 
chapter on Kitamura's method to compute the values of E, FS , 
Foy and FL/F, Chey characteristic functions). © LCHT2? isa 
WATFIV program (a listing and sample run follows), and is 
therefore initiated by 

SR *WATELY SCARDS = sLCFETI2.. 

As with ZAPP.DC, all results are printed out at the terminal. 
If one wanted to save the results ina file, one would 
follow the SCARDS command with SPRINT = an output file. 

The only input data required is the set of observations 
(half of one eclipse). The integer variable INORM, which 
is an input parameter, allows one to either normalize the 
Poh bscunve isingdeequation (32.5) CINORM = 1)Sor ce OMmiLesthis 
step (INORM # 1). The author's experience is that there is 
little difference in the values of E, Flu Bo and FL/F5 
when INORM is used. One must also enter the last two values 
of the phase 6 if INORM is to be used (in fact, these two 
values of 86 must be included regardless of whether or not 


INORM is used). Hence, the first input line reads: 


number of observations, 9 8 INORM 
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The observations in (0,2) pairs follow immediately. ULCrT2 
DEINESHOUt» aml ot the observations, followed by the values 
Of the transiorms (oF and Cady and finally the values of E, 
Foy Foy and FL/F 5: All computations are done using double 


precision (REAL*8) arithmetic. 
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Appendix 3 Programs for Kopal's Method 


There aré Six programs for Kopal"s method: TOTAL 
(total eclipses), ANNULAR (annular eclipses), PARTIAL (par- 
tial eclipses), EB.FS (least-squares Fourier analysis and 
moment determination), KAL (Kalman filter for moment deter- 
Mination), and ERROR (error analysis for orbital elements). 
Programs EB.FS and KAL have been adapted from the programs 
publishéedebysdunkeviens (1981) 2a Ali@six programs use the 
FORTRAN language, and are used in the following way: 

SHELTS= 150 

$R *FORTG SCARDS = PROGRAM NAME T = 1.0 

(computer responds with a message) 

$R - LOAD# 5 = INPUT DATA FILE 6 = -P 

(another message). 
Posesees ehewresultsp7elisty chestemporary L1ilee—-P. ) Lstingsnor 
all programs, as well as sample input and output files 
follow. The programs TOTAL, ANNULAR, and PARTIAL all re- 
Guire the moments of the light curve as input data. in all 
cases, a value of the limb darkening is required. ANNULAR 
and PARTIAL require the depths of both minima as input data. 
ANNULAR uses the iterative method devised by Jurkevich (1970, 
bg. 75) for determining k (See chapter on HS Herculis for a 
complete description. PARTIAL requires an initial estimate 
of k, which is best determined by forming (according to 
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AC Ne 
g, = : and g, = : 
Z AoA 4 AAG 
where G5 — 95 arc.) and J, =9, (arc). The parameters a and 
Cc, are defined as (xr, = radius of eclipsed star, r= radius 
Of eclipsing star): 
ee = eae cos il 
Bie Sie sy eal: 


One may determine a and om using the set of tables published 
by Demircan (1981, pp. 144-147). Another parameter b=l1-a 


may be used to find k: 
2E8. 
5 A 


This will be either greater than or less than one according 
to the eclipse type. To establish the eclipse type in the 
case of a partial eclipse, one may take advantage of a use- 
ful property of the Russell-Merrill y-functions 


vo° (n=0.8) > x -" (n=0-8) 


where 


(2 = value OL ev ab Minimum). 


The x-functions are easily computed with the equation 


esac Gay 
ilar Getic 
The Y=—tunctions, allow fon a quick determination of the 
eclipse type. Table A3.1 summarizes the input data required 


(in, the correct format) by TOTAL, ANNULAR, and PARTIAL. 
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Program EB.FS requires only the ascending branch of one 
eclipse as input data, along with some Out=OEr—eclipse 
points for the determination of Au. EB.FS requires the 
angle of external tangency, as well as a number of integer 
numbers for controlling the process of calculation. The 
format shown in the sample input file is the one which will 
be most commonly used. If the moments are not to be com- 
puted, IYN in dine 71 should be set equal to 2.) Thevoniy 
other parameters of interest are MI, the number of in- 
eclipse observations, MMAX, the number of out-of-eclipse 
observations, and NMAX, the number of Fourier coefficients 
to be used. These three numbers appear on line 1 in the 
following order: 

NMAX, MI, MMAX . 
The value of 6 at external tangency follows these numbers. 
Lines 3 need not be altered in most cases. The last two 
Mines of the EB-FS input file read: 

ONexteriads itangency.) e's nel pe. 

1.0EB-06, 
where 1 and 7 are digits controlling the number of moments 
to be calculated, and 1.0E-06 is the maximum allowable error 
im, the calculation of the moments. Program KAD is quite 
simple to use. The first line of the ipit wel Lencontains 
the value of q/r, which is denoted by T in the program (q/r= 
intensity of white noise/variance of variable), the value of 
1-2 at zero phase (estimated), and the number of observa- 
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Ty, 2-2 (0), munber. of observations. 
ine Gata ToOllow in (@,%) pairs, where © must be in radians. 
Program G€rror requires a ‘special command for execution. 
Since it incorporates in IMSL subroutine, the EOL owing 
command should be used: 

R= HOADH +7 IMSLULS 5 = aarout fille, 6 = =P. 


Tiew INGA eOLMat mequired aS. Showm, tm table A321, 


Lor 


mt ae 


aes a 


Le >). 
: - on 
mele «os 

pans > aio 


a) 

a= PH aaa = 2 
st. t& @ldew mie ; na r 
hee PL Pd bar ari 

a vs ica _ _ 


ru 


7A bi) Legh Fees 
af 


LG? 


lable Ao ee lnpuL t Ormatsefor Programs Using Kopal Method. 


TOTAL: (Neer aes A (format 4E15.6) 


Ghormmat, 2514147.) 


ANNUGAR: x, DOC, DIR," Is (format 4F14.7) 


A_,A 


5 91 Ay (format 4615.6) 


(NOter #DOCG=watte of © at minimum of occultation, 


DIR =avValue Ome at minimum of steransit):. 


PARTIAL X7DOC,aDIR, ki(estimated)) (format. 4F 14,7) 


Ags Ady Aye Ae (format 4615.6) 
eis OG 
ERROR: AAG , AAS 6 AAs, AA ys AAcr NUYPR (format. Sh l5.6, 13) 
ri, ro Na i,k, 1, , x (format SF1l4.7) 
FNS ay (format 3F14.7) 
fe) O O 
where: Eee = eransit 
Oc = Occul tation 
ib occultation 
NTYPE = 
anything else transit 
Liven eclipsed star 
r4 = eclipsing star 
X = re csc” i 
Y = bots 
@ = value of a at minimum of primary eclipse 
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PROGRAM TO ESTIMATE MOMENTS OF BINARY LIGHT CURVES FROM 170 


FOURIER COSINE SERIES APPROXIMATION OCF LIGHT LOSS. 
THIS PROGRAM ACAPTED (FROM ARTICLE BY JURKEVICH) BY 
De HOLMGREN FOR USE ON MTS FORTRANe NOVe 19'83¢ 


IMPLICIT REAL*8 (A—H,0—Z) 

DIMENSION TmC400) SEL (400s RUCI S615) s RIL G1ss1s)'sG Cl Se13)). 
~AER(13)9CS(13)sSN(13) 9X (13) »sAPHI (959) AMOM(13),AMOMER(13), 
=-ePs Cl3) 


DIMENSION OF 13 ON RU AND RL IMPLIES 13 TERMS IN THE REGRESSION 


EQUATION:10 HARMONICSsA09(=DELU),AND THE ABSOLUTE TERM. 


READ(Ss913) NMAXsMI +MMAX 
LS FORMAT CL 2.2213) 
READ(Ss14) THETAT 
14 FORMAT (F14.7) 
READ(5s499) 12,13 
499 FORMAT (2T2) 
THETATS 1S THE® POINT OF EXTERNAL CONTACT IN RADTANS. 
IF [I2=1 DO NOT COMPUTE INDIVIDUAL RESIDUALS FIRST TIME ARCUNDe. 
IF 13=1 (AUTOMATICALLY SETS 11=0) COMPUTE RESIDUALS. 
NMAX=HIGHEST HARMONIC USED» MI=NUMBER OF OBSERVATIONS WITHIN 
ECLIPSEs MMAX=NUMBER OF OBSERVATIONS OUTSIDE ECLIPSE. 


wWRITE(6,s800) 
800 FORMAT (1H1) 
WRITE(6s120) NMAXoMI » MMAX 
120 FORMAT (1X SHNMAX=12¢ 2X s3HMI=51 35 2X3 SHMMAX=I 3, 2X) 
WRETEG6s 121) 12513 
121 FORMAT (1X%5SH T2=12,5 2Xs3HI3=I12) 
WRITE(69122) THETAT 
122 FORMAT (1X5 7HTHETAT=F8-6//) 
NOBS=MI+MMAX 
DO 15 !=1,NOBS 
READ (S616) TCI). EL CT) 
16 FORMAT (2E15-26) 
1S CONTINUE 


DEFINE, CONSTANTS. 
PT=3-141592653589793D0 


I MAX=NMAX#2 
JMAX=NMA X+3 
REF ANG=PI/THETAT 


COMPUTE VEGSS9 or WwihGrir.s 
DO 53 I=1,MI 
EL(I)=1.-—-EL(1) 

35) (COIN MINUE 


COMPUTE SUM OF SQUARES GF LIGHT Loss 
De) 62 Soi 
TRY 1=TRY T+EL CT) *ECCI) 
se CONTINUE 


ZERO GUT RU MATRIX. 
DO 1 J=1sIMAX 
DO 2 If=1ls5J 
2 RUC JoI )=06 
RUC IJMAXs J) =06 
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TS CONTIUNUE 
é 
C COMPUTATION OF COSINES FOR K*TH OBSERVATION.’ 
11 DO 10 M=1,5MI 
ANG=REFANG*XT(M) 
SSN=DSINCANG ) 
CSN=OCOS(ANG ) 
SN(1)=SSN 
CS(1)=CSN 
DO 200 J=2sNMAX 
SNC JI=SN(CJ—1) * CSN#CS (J—-1) *¥SSN 
200 CS(JI)=CS(CJ—1) * CSN—SN(J—-1)*SSN 


C COMPUTATION OF INDIVIDUAL RESIDUALS IF NEEDED. 

LE Glh2— Py 206s 207.4 h02 

20E Y=0~6 
DO S01 I=2sNMAX 

SO1l Y=Y+Q( Is JMAX)*CS(I-1 ) 
Y=Y+02¢5*Q( 15 JIMAX)—Q( IMAX»JMAX) 
S1=EL(M)-Y 
S$3=S3+S1*S1 
WRITE (63128) T(M) sEL(M) sYsS1 

l28 FORMAT Cl0X%s0254165.5(5%D25.16)) 
I3=0 
GO EO 10 


C COMPUTATION QF COEFFICIENTS OF NORMAL EQUATIONS. 

207 XCU1)=0' 65 
NU=NMA X+1 
DO 201 J=2sNU 
X(JI=CS( J—1) 

201 CONTINUE 
X( JMAX )=—-EL(M) 
DO 203 J=1lsNU 
DO 202 I[=1,5J 

202 RUCJsI)=RUCJsT1)4+X(1)*X(4) 
RUC IJMAX» J) =RUC IMAX 9 J) +X (5S) ®X( MAX) 

203 CONTINUE 

LORE ONT NUE 

IFC l2sEQ@el)GGQ tO 41 
WRETECGs: P2795 3 

127 FORMAT (7Xs 3HE3=0232-16) 

41 DO SOO J=1sIMAX 
RUCIMAXs J )=—-2¢*RUC 51) 

S00 CONTINUE 
RUCIMAXs IMAX)=RUCIMAXs IMAX)+MMAX 
RUCIMAXs I MAX)=—2¢*RU(IMAX 91 ) 


cE SOLU ICN (ole NORMAL EQUATIONS. 

NP 1=JMAX 
CALL CRACOV(RU sNP1 sRL 9Q) 
WRITE(C E5124) 

124 FORMAT (/S/) 
WRITE G61 1h2)5)) 

125 FORMAT (2X%s1HJ 93Xs1lHI 312Xs2HRUs 24X%s 2@HRL 1 8Xs1LHI»s3Xel1lHJ, 

Xe liiG/ay 

DO SOS I=1,JMAX 
DO S04 J=I4.JMAX 
CURE ENC S) 5 TZ SIN A) 5 UIUC io ORIEL CIO 1) oI OA ONE tho J) 


123 (FORMAT Ci Xe l2s2X 91 2i03X%sD023616 3X eD250e1635X%s 120 2X%s [255% sD25e16) 
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seh 
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LSS! 
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LVS 
136 
rS7 
PSs 
LENS, 
140 
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SIG 
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iVisy2 
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CONTINUE 
WRITECEs506) 
FORMAT (7) 
CONTENUE 


CCOMPUTES THE SECOND. TERMAOPEQUATICN (15). 


c 


36 


TRY2=06 

DO 36 I=1lsIMAX 

TRY 2=TRY2+RL( JIMAX 4 1) *¥RLCUMAXoI ) 
CONTINUE 
DEN=MI+1—IMAX 
Sa=(URyl=TRY 2) 7 DEN 
ERRMEA=DSQRT(S2) 


C COMPUTE ERRORS OF UNKNOWNSe 
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DO 33 If=1,IMAX 

AER(I)=0-6 

CONTINUE 

DO 35 I=1,IMAX 

DO 32 J=IsIMAX 
AERCI)=AER(1)+Q( 1,45) *Q(I +d) 
AER(I)=ERRMEA*DSQRT(AER(I) ) 
CONTINUE 


C PRINT UNKNOWNS AND THETR ERRORS 
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WRITE(65124) 

WRITE( 6438) 

FORMAT (6X9 LHI 512Xs4HA( 1) s22Xs6HAER(I)DS) 
DO 37 IT=1sIMAX 

fA=it<—1 

WRITE (69 39)112Q(I5JMAX),AER( 1) 
FORMAT (S5X%e1T292(03X%sD023-16)) 
CONTINUE 

WRITE CEs 126)'S2 

FORMAT (7 7Xe3HS2=D23416//) 

DO 6186 L=1e2 

WRITE 6800) 

CONTINUE 


TEC §3—1)1025208s102 
T2=0 

S3=06 

WRITE(64110) 


110 FORMAT (20Xs4HT(M) 5 20XsSHEL(M) s22Xe1LHY¥s25Xs 2HS1 ) 
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READ(S 017) IYN 

FORMAT(1I1) 

IF CIYN«weEQe1) GO TO 612 

READ(S5 218) THETAsMUIN»sMUFIN,IINSIFIN 
FORMAT (F 14675412) 

READ (C5919) “EPS.Ct) 

FORMAT (E15¢6) 

DO 617 K=259 

EPS(CK)=EPS(CK—-1)*0.1D0 

CONTINUE 


C NOTE THAT MU AND I ARE OFFSETe TO RUN MU=0 CASE, 
C THETA IS IN RADIANS. 
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G 173 
OC 12 J=MUINsMUFIN 


TOL=EPS(J) 
00 12 I=TINsIFIN 
CG SELOW Jt 1S MUst? -tS THE RUNNING INDEX "IN “PHI T.sMU). 

J1=J-1 
I1=I-1 
WRITE C6s 1OSV I 151i, EPSCJ) 

103 FORMAT (1 Xs 3HMU=135 3X9 2HI=1393X »4HEPS=D236.16/) 
CALL PAT GhlsJ1s THETA sTGLs JJs PRI) 
WRITE C6. FO47 JT sills JJ s.PHr) 

1904 FORMAT (1X s3HMU=4 13 53X s2HI=135s3X 9 3HJJ=9 162 3Xs10HPHI (1 +MU)=D23016//) 
APHI(I»J)=PHI1 

27 GONTINUE 


C PRINT TABLE OF PHI*S. 
WRITE( 6800) 
WRITE( 65614) 
614 FORMAT (2X s2HMU s2X o LHI 9 SXs 1 OHAPHI (1 sMU) o11X s7HEPS(MU)/) 
DO 616 J=1,MUFIN 
DO,619 I=f. le IN 
Ji=J—1 
1Sy—1 
WRITE (62615) J1 oI1sAPHI(I1oJ)sEPS(J) 
615 FORMAT (2X%s7T2s2X%s1193X%sD2361693XsD8e2) 
619 CONTINUE 
WRITE(€,506) 
616 CONTINUE 


C COMPUTE MOMENTS FROM ASUBMU=SUM(CAPHI*NEUMANNNO *ASUBI ) 


IFIN=IMAX—1 

DO 300 J=1+MUFIN 

I=1 

SUM=APHI(I+J)*Q(1I+JMAX) 

DO 400 I[=2,IFIN 

SUM=SUM#2¢*APHI(1IsJ)*Q(1sJMAX) 
400 CONTINUE 

AMOM(J)=SUM 
300 CONTINUE 


C COMPUTE ERRORS IN AMOMCT ) 
DO 600 J=1,s,MUFIN 
SUM1=)0-e 
DO 602 I=1,IFIN 
SUM=0-6 
DO 601 K=1el 
601 SUM=SUM+APHI(K oJ) *Q( KI) 
SUM1=SUM1+SUMx*xSUM 
€02 CONTINUE 
AMOMER( J )=ERRMEA*DSQRT(SUM1 ) 
600 CONTINUE 


C PRINT A TABLE OF MOMENTS AND THEIR ERRORS 

WRITE( 65,800) 
WRITEC 6s 610) 

610 FORMAT (10X%s2HMUs1l12X%s 6HMOMENT >» 22X%s SHERRORZ ) 
DO €12 J=1,MUFIN 
MU=J-1 
WRITE (6s611)MUsAMOM( J) sAMOMERC J) 

E11 FORMAT (10X01234X%s023041634X%,0235016) 

62 CONTINUE 
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STOP 174 
END 

SUBROUTINE CRACCV( RU, NPl1 so RL>Q) 
IMPLICIT REAL*®8 (A-H,O-Z) 
DEMENS TONE RU CESB Sy) a Ria Gls 41S) sO 1S5t 3) 
LON PES Sie) ID 

N=NP1—1 

KL(NP1»NP1)=1.- 

RL( 141 )=DSQRTCRU(141)) 

DO SO K=1;3N 

KM1=K—-1 

DO 40 I=K,NP1 

SUM=06 

IF(KM1eNE20)GO TO 10 
RL(CIsKI=RUCIsKI/SRL(151) 

GO TO 40 

DO 30 J=1sKM1 
SUM=SUMF#FRL(I5J)*RL( Ks Jd) 
CONTINUE 

GEIGK es = Gierl GO Geo 
RL(TsKI=CRUCTsK)—SUM)/RL(K 5K) 
6@) Wo) 7o 
RL(IsK)=DSQRT(CRUCIsK)—SUM) 
CONTINUE 

CONTINUE 

CONTINUE 

DO 100 I[=1%5N 

OG elie ARIS Gl sale) 

CONTINUE 

Q(NP1.sNP1)=1le 

DOMES OM —iisiN 

K2=NP1-I 

DO 120 K=1.sK2 

SUM=06 

DO 110 J=1eK 

DI=J-1 

SUM=SUM#RL(I+KsI+01) *Q(1I,1+4+01) 
CONTINUE 

Q(1IsI+K)=—-SUM/RL( I+K,I+K) 
CONTINUE 

CONTINUE 

RETURN 

END 

SUBROUTINE PHI CI sMUs THETA, TOLs JJsPHI1) 
IMPLICIT REAL*8 (AH, U-Z) 
eA oiy, 

Pips l4iso2essasobDoO 

E=2 .71828182845D0 
ONE=FLOAT(-1) 

AMU=FLOAT(MU) 
FACT=( Ex *(1+MU))/PI 

CALL GAMMA(1e¢+AMUsGAM, IER) 
TSI=GAM/ (22e¢00**MU) 
AI=FLOAT(1) 

AIPI=AIx*PI 

TS So—CAMU +56 DIO) EZA Ae O 
TS4=(AMU412¢D90)/2-D0 

SUM=0-D0 
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PSeMUs ODOR CRE EVEN 


IODD=1 

Y=AMU 

Y=Y/2.-D0 
JY=SNGL(Y) 

Iya el XG) 
DIFF=Y—-FLOAT(IY) 


TEC DABS COIFF ) sbbsO061D=—10) LODD=0 


GENERAL, FORM OF SERIES FOLEOGWS « 


20 


J=J+1 
ITS=MU41—2*J 


FPCITS~.EQs0)GG TO 10 


AJ=FLCAT(J) 
TS2=20¢e *ASt1 eo 
ARG=TS2* THETA 
ARGSQ=ARG*ARG 


AS, 


TEMPNU= (ONE** (I +J) )*ARGSQ*DSIN(CARG)/(CARGSQ—-AIPTI*¥AIPI) 
IF ( (MU eGEe0) eANDe (Je GTe50))GO TO 21 
CALL GAMMA(TS3+4AJsGAM1,IER) 


-CALL GAMMACTS4—AJsGAM2,IER) 


FACT1=1¢/(GAM1 *GAM2) 


Go) ie Se 


21 AL=(CAMU+F1¢)/020¢*AJ) 
A2=(AMU~-10¢)/7(2e*AJ) 
ANUM=( 10e—-A1)**(J—-MU/2) 
DENA=(1¢4¢—A2) ** (Ut 1+MUZ2) 


DENB=AJ**( 1+MU) 


TEMP 1=ANUMZ (DENA*DENB ) 


AMPL=FACT*TEMPI1 


ANG=(AJ~AMUZ2¢+0¢500) *PI 


FACT1I=AMPL*DSIN( 


C COMMON SECTION. 


(S 


C PARAMETERS:XX — THE ARGUMENT OF THE GAMMA FUNCTION. 
GX — THE RESULTANT GAMMA FUNCTION VALUESe 


ia @ @ @@ 


ANG) 


30 TEMF=TEMPNU*FACTI 


10 


SUM=SUM+ TEMP 
IFC IOODeEQe1)GO 


TOy2Zo 


IF COABSCTEMP) -LE«TOL)IGO TO 10 


Ge) wer 2e@ 
CONTINUE 
JJ=J 
PHI1=TS1*SUM 
RETURN 

END 


SUBROUTINE GAMMA(CXX+GX>sIER) 
C SUBROUTINE TO COMPUTE VALUES OF GAMMA FUNCTION. 


eR => WOLERAINCE COIDEs 
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IMPLICIT REAL*8 
IF (XX eLEe34e5D0) 


NO TOLERANCE 
XX IS WITHIN 
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367 IF (X=2200)50950515 

268 10 IF(XeLEs2-eD0) GO TO 110 

seg 15 X=X-1.D0 

370 GX=GX* X 

371 GO TO 10 

372 S50 IF(X-1.00)6051205110 

373 C SEE IF X IS NEAR NEGATIVE INTEGER OR ZERO 
374 60 IF(XeGTeTOL)GO TO 80 

375 Kx 

376 Y=FLOAT(K)—X 

Bai IF(DABS(Y) eLEeTOL) GO TO 130 

378 IF((1eD0-Y)eLE-TOL) GO TO 130 

379 C X NOT NEAR NEGATIVE INTEGER OR ZERO. 

380 70 IF(XeGTeleD0) GO TO 110 

381 80 GX=GX/X__ 
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383 GO TO 70 

284 110 Y=X+1.D0 
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391 RETURN 
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PROGRAM TO ESTIMATE MOMENTS GF BINARIES BY 
KALMAN FILTER - EVEN MOMENTS AC2M) 6 

THIS 1S JURKEVICH*S PROGRAM, 

ADAPTED FOR USE GN MTS FORTRAN BY 


De 


3 
2 


4 


HOLMGRENs NOVe* 836 


IMPLICIT REAL*3 (A-H,O0-Z) 
DIMENSION F(400)..EL(400) 
INTEGER NOBS 

READ(Ss1) TeXsNOS8S 
FORMAT (2F 1427413) 

DO 2 T[=1s.NG8S 

READ(5,3) FCI) »sEL{ 1) 
FORMAT (2F 14.7) 

CONTINUE 

P=le 

P1l2=06 

P14=06 

P1i6=0-6. 

P18=0. 

X2=O06 

X4=O065 

X6=O06 

XB8=O06 

B=0-6 

DO 4 I[=1,»N08S 

Z=1--EL( 1) 

DIFZX=Z—-xX 

A=B 

B=F (1) 

SA=DSIN(A) 

SB=DSIN(B) 

C=B-A 

AM=P4+T*C 

SB2=S8*SB 

SA2=SA¥*SA 

SB4=SB2*SB2 

SA4=SA2*SA2 

SB6=SB2¥*SB4 

SB8=SB4*SB4 

DIF 2=(SB+SA) *( S8—-SA) 
SIG=SB2+SA2 

DIF4=SIG*eDIF2 

SSBSA=SB*SA 

OIF 6=D0IF2*( SIG-SBSA) *{ SIGt+tSBSA) 
DIF 8=(SB4+SA4) *SIG*DIF2 
AM12=P12+P*DIF24¢T¥*¥(C*SB24+S2(A)-S2(B)) 
AM14=P144P¥*¥DIF44*T*¥(C*SB44+S54(A)-S4(B)) 
AM16=P164+P*DIF64+T¥*¥(C*¥SB64+S6(A)—-S6(B) ) 
AM1 8=P18+P¥*DIF 8+T* (C*SB8+S8(A)—-S8{(B)) 
DM= be/(1 e+AM) 

P=AM*DM 

P12=AM12*DM 

P14=AM14*DM 

P16=AM16*0DM 

P18=AM18*DM 
X2=XK2+X¥*¥DLF24+P12*DIF ZX 
X4=K4+X*DITF44+P14*DOI1F 2X 
X6=X6+X*DLIF6O+P16*DIF ZX 
X8=X8+X*DIF8t+P1B*DIF ZX 
X=X+PxDOIFZX 
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61 WRITE(6s300)X22X45X6HxX8B 

[aye 300 FORMAT (1X, 3HAZ=E1105 s2XSZHA4S4=E1125+s2X3HAG=E11 05 52X3HAS=E1IL-S//) 
63 So 

64 END 

65 Ree Te WUENRE TE EGU) Sy22i0 08 )) 

66 Rue ve, Seo i 

ows S2=(X-SIN(X)*COS(X))/26 

68 RETURN 

69 END 

70 REAL FUNCTION S4(X) 

val REAL S4a9X9SXsCXx 

Y2 SX=SIN(X) 

eS CX=COS(X) 

74 S4=3e*® XS Be-( Be *¥SK/S/ Be t(SX¥SK¥ESXIS 4c) *CX 
75 RETURN 

76 END 

77 REAL FUNCTION S6{X) 

78 REAL S65X%sSX2TCX2535S5 

79 SX=SIN{(X) 

80 CX=COS(X) 

Bl S3=SKX*SXK*SX 

B2 SS5S=S3*SK*¥SX 

13) S6=SeFEXS16 0— (De ESKS LG 0 tHe KFSSBS24 0 ¢SS/SG6 ) ¥CX 
B84 RETURN 

‘NS, END 

86 REAL FUNCTION S8(X) 

87 REAL S8s Xe SXsCX35338S5%S7 

88 SX=SIN(X) 

89 CxX=cCOS{(xXx) 

99 S3=SX¥SK*SX 

91 S5=S53* SX*SX 

92 S7=S5*SK*SXK 

93 SB=35eFXS1L28 0 —( 350 FSKS/128 et 35 eo FSB/1LG2e Ht TAL*ESSS4 Be +tST/Be ) ¥CX 
94 RETURN 

95 END 
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PROGRAM TO USE K@PALY S FREQUENCY DOMAIN METHOD 185 
WOVCOMPUTE GEOMETRIC UEEEMENTS OF” AN@ECEI PSING eI NARY 

THiS V.ERSTPOND FOR. TOGA ECLIPSE S!=— PROGRAMM" TOTALS. 

De HOLMGRENs AUGe *8S6 

Rie WI SIS(D) VNC Ge (OR USE Wille Jesh oe Soc 


IMPLICIT REAL (A~Hs,G-Z) 
DI MENSTONY C.C3)% C803) 
READ IN DATA. 

READ(591)A0sA2 oA4 5 A6 

1 FORMAT (4E15-6) 
READ(Ss2)UsEL1 

2 FORMAT (2F14e7) 
CB(3)=A2/A0 
WRITEC 635) AODsA2sA4sA6 

S FORMAT (4E15e¢5) 
C(3)=CB( 3) 
CB(2)=(A0*A4—A2*%A2)/7( A0*ADN) 
CB(1L)=(CAOXAD¥A 6-3 0 FAD KAS EAS Fle FAZEALVKA2ZIS (AOXAOXAD) 
CBi1)=CB(1)/7CB( 2) 
C(2)=CBC2])*5 eX (3 0e-UIS (150-70 *U) 
CQLI=HACB( LIE CEC 1S e—-7Te KFUIS (30% (350-19 6%*U) ) 
D=(146-C(3) 2 *C(1)+Cl2) 
R1=C(1)*CC(1)70 
R1=SQRTCABS(R1)) 
R2=SQRT(ABS(C(2)7D)) 
SSi1=C(1)7D 
WRITE 65 3) 

SSE-GRMATH* ORBITAL ELEMENTS: § 73 
WRITE(E54) RlsR2sSSI 

4 FORMAT CPR1="eF 14071 Xs P*R2="oFl4e7oikXs*SINI¥*2=' Fl 4e 7) 
STGP 
END 
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TO COMPUTE GEOMETRIC ELEMENTS OF AN ECLIPSING BINARY 
THIS VERSION FOR ANNULAR ECLIPSES - PROGRAM *ANNULAR* « 
De HOLMGRENs AUGe *33% MODIFIED SEPT. 18°83 


(7) (ie) A 


REAL AOD sA2sA4,A659R14R255S515C(3)+,CB(3), 
—-Xe Y¥UsCC(2)e¢YYsINCL CTI sPI sTAUsDKeRKL 
—sRK9S1 sS2slL2e11 sALFOs ZsPIsXGeAMsFIeRK1.0GC,DITR 
P1=3.214159 
READ IN DATAs 
(ol CRs) Ure ints Geeks Spey Myer is 
DEPTEEOR Ie CCU TAR ION. Dil R tA sO sae 
TRANSIT. 
READ(Ss3) U,D00T,DTRsLil 
3 FORMAT (4F14-.7) 
READ(5 34) A02A2sA45AEB 
4&4 FORMAT(4E15-26) 
WRITE{(625) AQsA2s,A4sA6 
STE GRMAT CA0= 2 E1560 5 A295 £1 Ss Os PASH ee 15 eS “AD—" 561525) 
RK=065 
RKL=RK 
DO i 1=1,10 
TAU=206*(3e*ARSIN(SQRTORK ) I-03 0-4 e FRK IEC LOH Ce*RK) * 
—-SQRT(RK*(12—-RK)))/(3-*P1) 
VYVY=(3-*(12-U)+26.*U*¥TAU/S(RK¥RK))/4(3--U) 
RK=(12-DOTR)I/S(DOCKYY) 
RK=SQRT( RK) 
OK=ABS(RK—-RKL? 
IF (DK eLEe1-E-04)G0 TO 2 
1 RKL=RK 
2 CONTINUE 
RK=1e/RK 
G CONVERT FROM RUSSELL *S Ke 
C RK IS Ks K=RI/SR2S (K>1) 
X={( AD*¥ AD ¥AG6—3 2 FADFKAZ*XA4G4+2 oe KALKACKALIS(ADKXA4—-A2Z2¥*A2Q) 
F2=(1542-76*U)/(52*(3-2-U) ) 
F4=34%*§ (350-194 *U) S( 356 *(3--U)) 
X=X*(AO¥ (1Le2tF2*RK*¥RK)D-L1) 
X=XSCAO*FAOK (Let be XF2AERKERK FF4O%* (REEKG4 62) J—-3e FADKL I KC LEFF AOXRK ERK) 
—+t2e*LI¥*LI1) 
Y=(X*L1-A2)/7A0 
C X=R2*¥*¥2 *¥CSCI**2, Y=COT [**2.6 
WRITE(69183) RK oeXeVYoelL ls YY 
18 FORMAT (* K=" .Fl4e7s*X=* sF1l4e7s* Y=" sFl4e7e*Li=* ,Fl4ae7 
—s'*VYY=",F14e/7) 
C FIND R1eR2e1INCLINATIONe 
IF(X el Te 04) X=ABS(X) 
IF (Y eLEs OoVINCL=1.25708 ot STOP 
[EV Pellet Os) GO TO) 16 Ge END 
Y=SQRT(Y) SIslo) (Oe MS jus 
INCL=ATAN(1-/7Y) 
16 WRITE(6, 30) 
30 FORMAT (*Y¥<=O54 [=904 ASSUMED .? ) 
R2=X*SINCINCLI*SIN( INCL) 
R2=SQRT(R2) 
R1=R2*RK 
WRITE( 6921) RI »R25TNCL 
21 FORMAT (®R1=* sF14e7 oP R2=" oF 1407s" INCL=* oF 14-07) 
WRITE(6,22) 
22 FORMAT(* THESE ARE FINAL ELEMENTS ® ) 
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GRAM TO” USE KOPAL*S FREQUENCY DOMAIN METHOD 

COVEPUPE GEOMETR TG SEEaMeENTS SOR VANS EGET PS ENG wet NAIRY 

S VERSTON FOR PARTIAL SCE ILPSES — PROGRAM “PARTI AL * .« 
HCLMGRENs MARe4!* 836 


IMPLICIT REAL CAH ,0-Z) 

REAL INGE ser 

PI=3e14159 
D IN DATA. 

READ(553) UsBGCsDTRIRK 

FORMAT (4F14-7) 

READ(5s4) A0sA2sA4sA6 

FORMAT (4E15¢6) 

WRITE( 695) AQsA2sA4s3A6 

FORMAT ( *AD=* ELS e691 Xs *AZ=*sE1LS5e691X 9 *A4G=* E15 e651Xs* AG 
DO 100 I=l1>s2 

RK2=RK ¥RK 

IFC RK eLT ele DALFAH=12=D0C+(1e—DTR)ISRKeS 

IF CRK eGE ele ) ALFA=1 e--DTRH( 126-DGC)*RK2 

L1I=AO/ ALFA 

IS Ks K=R1/R2 

X= (A0* AD *X¥A6—3 6 KAD XAZKAGH2CKEALEAC*EALZIS( AD*X¥A4G—AL¥*AC) 
F2=(15 0«-7e*U)/(5e¢*(30e-U) ) 

F4=3 0¥ (35 0-196 *UI/(35-%( 2e-U)) 
X=X*(AO*¥(1LetF2*RK*XRK )—L1 ) 


Rs je7) 


U6 fe oat?) 


X=XS(CAO*KAO¥F( Le t3be*X¥FA*X*RK*XRK FF 4*( RK¥*4 0) )—3e X¥AOXLI*( Let 2*RK*RK ) 


—-+t2exlLl *iLil ) 
Y=(X*L1—-A2)7A0 

2*xe *CSCI**2, Y=COT 1**2-6 
WRITE(6918) RKoXsYoLll 


FORMAT C8 K=* 6 Fl 407 1X oT X=8 Fl 4 eo Sol Xo* YH" oFi de 7s 1X s*Ll=" sFl4e7) 


D R1sR2sINCLINATIUCNe’ 

IF(X el Te Of VX=ABS(X) 

IFC Y elEes O«VINCL=1-5708 

ERGY velGce O's.) “GOV EO CLS 

YeSGRT CY.) 

INCL=ATAN(1¢/Y) 

CONTINUE 

R2=X*SINCINCL) *SINCINCL) 

R2=SQRT(R2) 

R1=R2*RK 

WRITE( 6921) R1lsR2sINCL 

FORMAT (#®R1=* gFl4e 7 s*® RE=* sFlI4 e075 * INCL ="5F14e27) 
RK=(A4—-AO*Y*¥Y)/S (CLI ¥*F O*X¥X DFC 20 *(Y/X Ile d)/F 2 
RK=SQRT(ABS(RK)) 

CONTINUc 

WRIMECGs 22) 

FGRMAT CATHESE ARE FINAL ELEMENTS * 2 

Sia 

END 


= © 


Amadeo 


a. 
* IN TTRAS, ‘ 


© o> i 


ML « Cilia *SHn* + Rle ae PsA ea « “2Ga" 174hMROR z 


ChA* SA~P AOA INES mo neygue os 68 Oe 4 DAM ak SAGA SOE DSH mi 


‘er umine.i~( (, bean bre Foal englee ae. b+. i) eOReOR eX 


ae YO bel ge doit ah ile ely (ae PTAMRON Gt : 


$ew8dSe) 3 
ae ‘saad 0a t4eaiatenw 


. Se4=1 064 80 

AUR MSHA 

aw DA dae se Jd with ive Tie WEP 
SAR YD DOTY el Aa deo ie Se BATE 
kA OWE 

GANDAes «& 21 RO 


C(0= Skt BN tube epee | 
( (ue ty, SE PAGS etm. SES ESR 
Pit 2 Ota S Se ad USAT ORSA ~—s 


4.2" 11%eSt 

PSV ( heed JOA jee 
.24¢8, Tibe’ .-S®ensFow Hewsrnak 3 
Ly Vga Fi (O43 03TI AG a 


ove nth *L DADs Ae dH OVA O 
(SBP AE «Ad «tds Ce te 
atest « eit iaet w Bde vba as 
Si OY Ge FeO ode VPA 


iv) Tesexy ye 
Ne pee ae ANG 
suecrewi> oi 
foam tae (Jor ke Lae xe ' 
(cH 7 Tave=esa 


| xae Sonn 
Lin SO. Cones PROM 
(Ne Mian, *SID PA Ue O14. (See ppm ape vs 7 1S 


SANK ed UN NY ) ee bo KES ee pene 


ad INLD eS 


mi 


OMNoOUF UNE 


10 


@) @ Oe oo © 


(a) 


C Seu 


SE 


PROGRAM TO COMPUTE UNCERTAINTIES IN 188 
PUEMENTS (GIVEN ERROGS OF MOMENTS AND VALUES oF 
EUEMENTSs (USES IMSL ROUTINE LEGQTIF TO SOLVE 

A SYSTEM OF EQUATIONS. 

BY De HOLMGREN. FEB. 64. 


10 RUNS SR =-LOAD#SE*INSLLIB S=<INPUT Frile> c=—P. 


IMPLICIT PREAU CA=H's 0-2) 
REAL INCU 
DIMENSION A393) +8(3)9WKAREA(10) 


FREAD IN MOMENTSsELEMENTSs AND ERKCRSo 


UP INITIAL ESTIMATE CF DELTA K**¥26 
READ(Ss1) DAOTRs DAD0OCs DA? sDA4S sDAE SNTYPE 
FORMAT (S5E15-6313) 
READ(5Ss2)R1sR29X%sYsINClLosAKsELl sU 
FORMAT (8F14.7) 

READ(Se15) AQOCsAONTRsALFO 

FORMAT (3Fi4e7) 

AK2=AK * AK 

DK2=1-E-03 

IF CNTYPE=-1)39403 
DALF=DAOTREAK2¥*¥DANOCX* (DAGOCSANOC HD K2/AK2) 
DL1=(DAOOC/ANDOC—DALF/ALFO) *EL1 

GC ye S 
DALF=DA0GC+(DADTR/AK 2) *(DAOTR/SADTR#DK2SAK 2 ) 
DL1I=(DAOTR/SAOTR=DALF/SALFO) *EL1 
CONTINUE 

UP EQUATIONS FOR DXsDYsDK**2e 

IF (CNYTPE-1)657 96 

AO=AO0OTR 

DAQ=DAOTR 

GG} 1G) 

AO=AO00C 

DAQ=DA0CCT 

CONT INUVE 

A(1.1)=EL1 

A(1ls2)=-A0 

Al1ls3)=0-64 

B( i J=DALN?]+Y*ODAU—X¥*DLI1 
F2=(150e-7e*U)S(Se*(30e~U) ) 

F423 0% (350-19 e XU) S (S50 *( Ser) ) 

A( 291 )=20e*ELI*( XH YtF2*AK2*X) 
Al2s2) =—-2e* (ELI *X—AD*¥Y ) 
A(2s3)=ELIL*F2Q* X* X 

B(2)=DA4-Y*Y*DAO 

B(2Z)=B( 2) ~ (CKKX— De KX*¥YHEFS]KAK2*X¥*¥X) X¥DOLI 
XMY=X~ Y 

FK2Z=F2*AKe 

FKG=F4*AKEx*AK2 

X2=X*¥X 

X3=X*X2 

Y2=Y*Y 

X Y= X*XY 
A(3e1)=XMY*¥XMY FF Ke ¥X *¥( Ge K¥XKM VN eK YI HFK OX 2 
A( 301 )=3-*ELI1*A(351) 

A(3s2) =EL1* ( X2—2 6 *XYtFKE*¥X2)+A0*Y2 

A( 332) =—-34-*A(332) 

A( 393) =3 eo ®F2*¥K2*XMYH 2e X¥F4¥* AK 2% XB 

A( 3,3) =EL1*A(353) 
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B(3)=DA6t+Y 2*Y*DAO 

Z=X*(CX2—-3e¢*¥XVYHIB*X*Y2) 

Z=Ltge *¥FKA*X2KX KMYFtFK4 *XS 
BC3IJ=8(3)—Z*DL1 

M=1 

N=3 

ITA=N 

IDOGT=0 

CALL LEUTIFCAsMsNoIAsBsIDGTsWwKARCAS IER) 


CesuaVe SCuUATLTONS ce X RETURNED SiN Bie 
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10 FORMAT (*DOCOS I1*®¥ *2=",E15 e690 1X5" DR2*X*2='5 ELS ee S91X%s*ORI**2=! 


WRETE C659) —CB CI )rs) f=1.4.3) 


FORMAT (#§DX="*sElSeGolXs*DY="sELSeEslXsa* DKXKXO=% FELSe ES) 


C53; —COSGINCE) COS CINCE) 
Ssil=le-—csl 

OX=6(1) 

DY=8(2) 

DK2=B( 3) 

DCSI=CY*CSI*SSI/SY 
DR2Z=R2*R2*(DX/X=-OCSI/ZSSI) 
DRY=RI¥*¥RI*OCDK2/S AK24tDRES(CR2*R2) ) 
WRITE(6s10)0CSI +DR2sDR1 


C RE-EVALUATE Dtl 


11 


l2 


is 


14 


TPONTY PEL) Tis 12.i1 
DALF=DAOTR+t+AKe2 *DAU OC (DAQOC/S AN OC +DK 27 AKe ) 
DLL=(DAOOC/AO0DOC=DALF/SALFO) ¥EL1 

Go wa Wes 
DALF=DAO0OC+( DAOTRK/ZAK2 ) *(DAOTK/SAOTR@=DK2/7AK2 ) 
DL1=(CDAOTR/ZAOTR-DALF/ALFO)*¥ELI 

CONTINUE 

WRITE(6514)0L1 

FORMAT (* NEW OL1="5E15¢6) 

STGP 

END 
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1 0.9020E 00504648255E-O01 6 0e37023E-025 06 34278E-03 
2 02630-90925 
ENDOOre ra Le 


input Eile, total eclipse (W Delphini), 


i 026200902030 47050202-9020 
2B 023338 02 01874540.0015258,0.2.0001545 
END) OF File 


Input file, annular eclipse (HS Herculis). 
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Appendix 4 The WINK8 Program 


The light curve synthesis program WINK8 is a modified 
version of Wood's original WINK program published in 1972 
(see Wood (1972)). WINK8 incorporates a number of improve- 
ments, only two of which will be mentioned here. The ics 
of these is the use of model stellar atmospheres, since 
these are more realistic than the blackbody approximation. 
In the program, these model atmospheres are used as ckigexehe 
data and are stored as tables with temperature, surface 
gravity, and wavelength as the arguments. The tabular quan— 
ity iS MesnOLlalkesubeace GlUx. nesses iluxes “are used to 
compute the intensity at each point on the star. An inter- 
polation, Youtine 1s used to derive the flux, which is then 
converted to an intensity (assuming the limb darkening to 
be a known quantity). The second major improvement incor- 
porated in WINK8 is the way in which the reflection effect 
is computed. Rather than using an explicit numerical imnte= 
gration scheme, WINK8 uses a series of approximations to 
compute certain quantities required in the calculation of 
the reflection effect. 

WINK8 is not difficult to use. Due to its great size 
(1857 Lines), 2 fs not reproduced here, but a sample data 
file is shown. The data file is composed of four distinct 
Darts. Thepesustecevencecn. lines Gontain 1npuc parameters 
Belabing toute sOroleeand eto the sstars themselves. The 


first line must contain the word 'WINK' in columns 1-4, and 
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any other alphabetic data in the remaining columns. Line 2 
must be blank: <bLines 3—-Lovcontain* the orbital and ‘stellar 
parameters. There are 19 parameters which may be treated as 
variables, but not all of these should be treated as vari- 
ables at the same time. In most cases to be encountered, 
the orbital and stellar parameters to be given specific 
values are those shown in lines 3-16. The important para- 
meters and their descriptions (along with any necessary 
comments) may be found in table A4.1. 

The “-1 0" must follow the last input parameter. Each 
input parameter has an index I which appears in columns 1 
and 2. MLE input parameters 1-9 vare not Known accurately, 
one need not specify any of them since WINK8 has a set of 
default parameters. However, such things as the wavelength 
Of Obsenpvatilon should) be specttteds = Lines 18-30 contain the 
model atmosphere data. The format shown is the one which 
must, bemused tinal cases] lines 16-24 Napply® tomtie primary 
star, and lines 25-31 to the secondary. The logarithm of 
the suntace gravity appeacs in lines U3 and 25,7 and at the 
Cndeo tel anece lOpe ie Joe aco eo OMAR G see ea Een 
peratures appear in lines 19, 22, 26 and 29. These tempera- 
tures are divided into two ranges: low (4000°K- 9500°K) and 
high (11000°K— 40000°K). These temperatures are followed by 
two lines of fluxes which correspond to the given tempera- 
tures. The wavelengths bracketing the observation wave- 
length are the second to last numbers on the lines CoOncadmLiag 


the fluxes. The model atmosphere data must be followed by 
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a blank line. All succeeding lines, except the last two in 
thee eile, contain the inpue data in the form: 

phase, magnitude or luminosity, observation weight. 
The phase may be in either radians or revolutions. The "=i" 
in the second last line must be present (anywhere in the 
first 10 columns). The last line is known as the "T/F card". 
This line tells WINK8 which of the 19 variable parameters 
are to be treated as variables and which are to be treated 
as constants. A "1" denotes that the parameter is to be 
varied, while a "0" denotes a constant. Note that each num- 
ber is separated by a space. 

To facilitate the construction of a data file, a 
demonstration file, known as DWINK, may be copied to an 
empty file and then altered to suit the star being analyzed. 

WINK8 will produce several pages of output, so an 
example is not shown here. However, the outcpue 2s entirely 
self-explanatory, so the user should not have any eae ybled ks 


in understanding it. 
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Table wie. l Important WINKS Input Parameters 
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digg : magnitude at quadrature equal to. 00) 2a 
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